program coordinates C C Version 3 - April 2007 - incorporates fits for LRIS-R and LRIS-B C of early 2007 with the ADC in and with the ADC out. Incorporates C a check of position angle as determined by the reerence star C X,Y(CCD) positions if there are 10 or more reference stars. C If the deduced angle differs from the position angle from C the telescope control system (keyword ROTPOSN = 90.0), there C is an option to redo the calculations using the new PA. C C Version 3.1 - Feb 2008, fixed quadrant problem with PA, C added 3 line section to subroutine solve_matrix.f C C July 2007: Compile: f77 -o coordinates coordinates.f C C C This program enables you to take (x,y) values (in units C of pixels on the LRIS CCD) for a list of stars and turn them into C RA and Dec on the sky. If you have good astrometry, a set of C reference stars of known RA and Dec can be used, otherwise the C coordinates of the center of the field as given in the image C header are used. Either Figaro or Fits images can be used. C Input is in the form of a file which contains the image name, C and x,y,ra,dec,epoch for the reference stars, and various C parameters for the object stars. The header of the image file is C read by the program to extract various other necessary parameters, C such as the hour angle, elevation, etc. C C There are two output files, coordinates.list, which is a file C suitable for use with Autoslit, the LRIS slit-mask assignment C program, and arcsec.list, which gives offsets N and E of the C field center in arc-sec for the objects before the correction C for atmospheric refraction is applied. C C A more complete description of the program can C be found in the help file, coordinates.doc. C C To use in the stand-alone mode independent of Figaro, change C "subroutine" on the first line to "program", i.e. "program coordinates". C This must be compiled as follows: (Note the -xl flag !) C C f77 -xl -o coordinates coordinates.f fake_dta_fig.f C C The "fake_dta_fig.f" file contains dummy routines to replace the C Figaro dta library calls. It must be linked in if "coordinates" C is to be run without Figaro. C C The -xl flag is part of the standard Sun Figaro compilation scripts. C There are no Figaro parameters used for this routine. If using with C Figaro, copy any routine's connection file, edit out all but one, C and use that to create a par file. If you don't do this, you will C get an error message out of the PAR subroutine library about C a missing file when you start to run it, which you can ignore. C C The program uses routines from Pat Wallace's SLALIB library. C For convenience, the source code for the routines required C is included in this file. C C The character string parsing routines, except those in SLALIB, C are from Dru Phillips, who wrote the input user interface for C Fits files. The Figaro header code was added by Judy Cohen, C who wrote the remainder of the program. C C J. Cohen - Caltech - Oct. 1994 C C WARNING: For simplicity, this version does not precess the C reference stars to the same epoch. It assumes that the epoch C of the coordinates of all the reference stars is the same as C the epoch of the coordinates of the first reference star. C C Modifications: C Correction of a bug found by A.Phillips of Lick. angle becomes C -angle in the rotation section. (April 18, 1995) C C Sep. 1996 C Added parameter "fit_flag" (an integer) C This parameter takes care of the difference in scale between C the original LRIS dewar and the Dewar installed in July 1996 C with the HIRES CCD and the meniscous field flattener designed C by Jim McCarthy. C Nov. 1996 - added full astrometric fits for the new LRIS Dewar C to subroutine transform. C July 1997 - fixed bug in section that write the output file. The C dec sign was incorrect for dec between 0(-epsilon) and -1 deg. C Feb 1998 - fixed decoding of date, reversed month and day C Nov. 2001 - C Decoding of date is wrong again. C This seems to have changed with year 2000, and also with C LRIS blue added, so instead I type it in. C C April to Aug. 2007 - C option -xl no longer supported by Sun fortran compiler, I removed it C from the compilation line suggested in earlier versions of C COORDINATES and in the documentation (doc) file. In this case C the record length must be set to 80. C fixed up DATE read C added check of PA between ROTPOSN from instrument rotator and C that determined from reference stars, if 10 or more ref. stars C are given in the input file C added new fits and split subroutine transform_blue and transform_red. C FIGARO files section has not been updated C C Note re LINUX (July 2007) C to compile in LINUX, modify the line C if (vstring(1:1).eq. '\\') then C to C if (vstring(1:1).eq. '\') then C C C implicit none C C Common array for x,y(CCD pixels into arcsec) and x,y RA,Dec into arcsec C from center of field, for determining true PA C integer rec_len, max_ref, max_obj, infile, image_lu C....NB: SET rec_len=80 (standard) or 20 (VMS or -xl switch) !!! had to remove -xl flag * parameter (rec_len=20) parameter (rec_len=80) parameter (max_ref=22, max_obj=1000) parameter (infile=11, image_lu=12) C....DEFAULTS: real def_epoch real def_mag integer def_pri parameter (def_epoch=2000., def_mag=20., def_pri=999) real*8 pi PARAMETER (PI = 3.14159265358979324) C....This are the parsing functions, and the variables used by them real get_real real*8 get_time integer get_int integer d_parse integer i1, i2, narg, len character*80 line character*6 lnfmt character*80 vstring C....These are the variables which contain the desired info C.. epoch = epoch for output coords C.. ra0, dec0, epoch0 = coords of LRIS image C.. rotang = rotation angle of LRIS image C.. refra, refdec, refepoch = coords of astrom. refs (arrays) C.. refx, refy = CCD coords of astrom. refs (arrays) C.. x, y = CCD coords of selected objects (arrays) C.. mag, priority, id = optional info for selected objects (arrays) C.....nref, nobj = number of ref, selected objects real angle,xarc(max_obj),yarc(max_obj),equinox real xarc_ref(max_ref),yarc_ref(max_ref) real north_arc_obj(max_obj), east_arc_obj(max_obj) real north_arc_ref(max_ref), east_arc_ref(max_obj) real epoch0, rotang, prepix, preline,date_install real blue_install ! date of installation of LRIS-B, added July 2007 real*8 ra0, dec0, ara_obj(max_obj), adec_obj(max_obj) real*8 ara_ref(max_ref),adec_ref(max_ref), rpas real*8 true_ra_center, true_dec_center real*8 true_ra_obj(max_obj), true_dec_obj(max_obj) real*8 true_ra_ref(max_ref), true_dec_ref(max_ref) real*8 refra(max_ref), refdec(max_ref),delta_ra(max_obj) real*8 elevation, hour_angle, delta_dec(max_obj) real*8 ra_average,dec_average,sigma_ra,sigma_dec real*8 dra_max, dra_min, ddec_max, ddec_min,dsum1, dsum2 real*8 center_obj_ra, center_obj_dec, dvalue, dvalue1 real refx(max_ref), refy(max_ref), refepoch(max_ref) real x(max_obj), y(max_obj) real mag(max_obj), date, value integer priority(max_obj), iarray_ra(4,max_obj),icenter_ra(4) integer icenter_dec(4), imonth, iday, iyear integer iarray_dec(4,max_obj), int_array(4),ichoice, ival integer ref_ok(max_ref),fit_flag integer icheck, iblue, iadc ! added July 2007 integer ione,itwo,itot, namp ! added July 2007 for binning check character*10 id(max_obj) character*80 string, string1 character*20 string20 character*64 zimage, image character*3 string3 character*2 string2 character*1 cvalue_ra(max_obj),cvalue_dec(max_obj) character*1 cvalue_center_ra, cvalue_center_dec,dchar integer nref, nobj, istatus, nstart, nref_bad C C Variables added for check of position angle (code from C program /scr/jlc/keckastrom/matrix.f C real*8 xtrue(max_obj),ytrue(max_obj),xxx(max_obj) real*8 yyy(max_obj) real deltax, deltay, theta_image, pa_tol C....These are some work variables character*8 keyword integer i, j integer in_stat integer section logical exists integer BACKSLASH C....Define the format; should be written with line-length given in parameters data lnfmt /'(A080)'/ C....Define the backslash character as ASCII code for VMS-style fortran data BACKSLASH /92/ C....Miscellaneous setup: zero counters, clear id-strings nobj = 0 nref = 0 do i = 1, max_obj id(i) = ' ' end do do i=1,max_ref ref_ok(i)=0 end do iblue = -99 ! LRIS-B or LRIS-R blue_install = 2001.5 ! Guess, based on LRIS News from Keck Web Site pa_tol = 0.15 ! tolerance for position angle solved vs that of ROTPOSN C C which version C write(*,*)' ' write(*,*)' This is version 3.0 of July 2007 ' write(*,*)' Coordinates, version 3.0, July 2007, includes ADC' write(*,*)' old version 2.0, Nov. 1996, modified ' write(*,*)' to handle the new and old LRIS Dewars.' write(*,*)' ' C C Default, ADC is on C iadc=1 C....Open the input file write(*,*) : ' What is the name of the input file with objects etc. ?' read(*,150)string 150 format(a50) open (infile, file=string, + status='old', access='sequential', err=910) C....read and parse the first line with the name of the LRIS image read (infile, lnfmt) line i1 = 1 i2 = 80 narg = d_parse (line, i1, i2, vstring, len) C....Open the image inquire (file=vstring(1:len), exist=exists) if (.not.exists) go to 916 open (image_lu, file=vstring(1:len), recl=rec_len, + status='old', access='direct', form='unformatted', err=911) C....Read in first line; check for FITS file read (image_lu, rec=1, iostat=in_stat, err=913) line(1:80) i1 = 11 i2 = 80 keyword = line(1:8) * write(*,993) keyword ! DEBUG * 993 format(a8) * write(*,994) line ! DEBUG * 994 format(a80) if (keyword.ne.'SIMPLE ') then close(image_lu) go to 200 end if i = 1 write(*,*)' this is a FITS file ' ! DEBUG * write(*,150)string write(*,*)' ' C....Read in next line; parse and check for relevant keywords 100 i = i + 1 read (image_lu, rec=i, iostat=in_stat, err=913) line(1:80) i1 = 11 i2 = 80 ! was 80 before keyword = line(1:8) * write(*,995)line ! debug 995 format(a80) if (keyword.eq.'END ') then icheck=0 if (iblue.eq.0) then write(*,*)' LRIS-R in use ' if (string20(1:9).eq.'TWOAMPTOP') then namp=2 write(*,*)namp, : ' number of readout amps with LRIS-R' else write(*,*)'PROBLEM - not standard TWOAMPTOP readout' write(*,*)' input number of amplifiers used' read(*,*)namp end if icheck = 1 end if if (iblue.eq.1) then write(*,*)' LRIS-B in use ' namp=4 icheck=1 if (date.lt. blue_install) then write(*,*)date, blue_install write(*,*)' PROBLEM - no LRIS-B fit available ' write(*,*)' prior to upgrade to 2 CCD mosaic' stop end if write(*,*)namp,prepix,' Assumed namp, prepix' end if if (icheck.eq.0) then write(*,*)' LRIS Red or Blue not correctly defined' stop end if write(*,*)' ' write(*,*)' ' write(*,*)' end of reading FITS header ' write(*,*)' ' write(*,*)' ' write(*,*)' ' go to 190 end if C C For parsing the RA, DEC, and HOUR ANGLE, which are stored as C character strings, we use some Keck specific code. First, for Fits C file cards, we remove the delineators, ":" and "'". Then we use C Pat Wallace's Slalib library routine DAFIN, which converts a C free-format character string from deg, arcmin, sec to radians. C C This is fine for DEC. For RA we then need to multiply by 15.0 C For HOUR ANGLE, we also need to multiply by 15.0 C C Phillips' original code uses: C dec0 = get_time (line, i1, i2, narg, 0.) C but I prefer the SLALIB library routines. They are better tested C and do not appear to have the compiler dependences that the original C versions of Phillips' character string parsing routines had. C C Note that consecutive calls to DAFIN working on different strings C must use different values of the initial search parameter, otherwise C the string is searched for a second angle further into the string. C C We use iblue = 1 for LRIS-Blue, iblue=0 for LRIS-Red C and distinguish them using the keyword INSTRUME= 'LRIS ' for LRIS-R C and INSTRUME= 'LRISBLUE' C C At the present time the ADC configuration cannot be determined C from the FITS header keywords. We read it in. C C if (keyword.eq. 'INSTRUME') then string20(1:20)=line(12:31) write(*,995)line(1:80) * write(*,999)string20(1:20) if (string20(1:6).eq. 'LRIS ') then iblue=0 else iblue=1 end if write(*,*)' ' end if if (keyword.eq.'RA ') then do j=1,10 line(j:j)=' ' end do do j=11,80 if (line(j:j) .eq. ':') line(j:j)=' ' if (line(j:j) .eq. '''') line(j:j)=' ' end do * write(*,995)line(1:80),' RA ' ! debug nstart=2 call sla_dafin(line,nstart,ra0,istatus) write(*,*)ra0,' RA, decoded ' ! debug ra0=ra0*15.0d0 write(*,*)' ' go to 100 end if if (keyword.eq.'DEC ') then do j=10,80 if (line(j:j) .eq. ':') line(j:j)=' ' if (line(j:j) .eq. '''') line(j:j)=' ' end do * write(*,995)line(1:80),' Dec ' ! debug nstart=10 call sla_dafin(line,nstart,dec0,istatus) write(*,*)dec0,' Dec, decoded ' ! debug write(*,*)' ' go to 100 end if if (keyword.eq.'HA ') then do j=11,80 if (line(j:j) .eq. ':') line(j:j)=' ' if (line(j:j) .eq. '''') line(j:j)=' ' end do * write(*,995)line(1:80),' HA ' ! debug nstart=11 call sla_dafin(line,nstart,hour_angle,istatus) write(*,*)hour_angle,' HA, decoded ' ! debug hour_angle=hour_angle*15.0d0 write(*,*)' ' go to 100 end if C C Read the amp configuration, for LRIS-R, 2 amp read, C IMTYPE = 'TWOAMPTOP' (LRIS-R) C IMTYPE = 'ONEAMP ' (LRIS-B) (old) C IMTYPE = 'n MOSAIC' (LRIS-B) (now) C AMPLIST = '2,1,0,0 ' (on 2007 images, not on older images) C BINNING = '1, 1 ' C if (keyword.eq.'BINNING ') then line(1:11)=' ' do j=9,50 if (line(j:j) .eq. ',') line(j:j)=' ' end do * write(*,995)line(1:80),' Binning ' ! debug read(line,*)ione,itwo write(*,*)ione,itwo,' Binning ' write(*,*)' ' itot=0 if (ione.eq.1) itot=itot+1 if (itwo.eq.1) itot=itot+1 if (itot.lt.2) then write(*,*)' Binning is not 1x1 ' stop end if end if if (keyword.eq.'IMTYPE ') then write(*,995)line(1:80) line(1:11)=' ' string20(1:20)=line(12:31) * write(*,999)string20(1:20) write(*,*)' ' end if 999 format(a20) C C The Fits format for date was C DATE = '15-09-94' C It appears that in 2000 the DATE format was changed C and it is now DATE = '2001-05-16' C So we have a complicated section to get the date correct. C Note that this is the true equinox of the coordinates in this C LRIS image. C if (keyword.eq.'DATE ') then do j=11,80 if (line(j:j) .eq. '''') line(j:j)=' ' if (line(j:j) .eq. '-') line(j:j)=' ' end do * write(*,995)line(1:80),' Date ' ! debug nstart=11 call sla_intin(line,nstart,iday,istatus) if (istatus.ne.0) write(*,*)' Error reading day ' call sla_intin(line,nstart,imonth,istatus) if (istatus.ne.0) write(*,*)' Error reading month ' call sla_intin(line,nstart,iyear,istatus) if (istatus.ne.0) write(*,*)' Error reading year ' if (iyear.gt.90 .and. iyear.lt.2000) then date=iyear+1900.0+float(imonth)/12.+float(iday)/365. else date = float(iday) + float(imonth)/12.0 + float(iyear)/365.0 end if equinox=date write(*,991)date write(*,*)' ' end if 991 format(f13.6,' date ') if (keyword.eq.'EQUINOX ') then epoch0 = get_real (line, i1, i2, narg, 0.) * write(*,995)line(1:80) write(*,*)epoch0,' Equinox ' write(*,*)' ' end if if (keyword.eq.'EL')then elevation = get_real (line, i1, i2, narg, 0.) * write(*,995)line(1:80) ! debug write(*,*)elevation,' Elevation ' write(*,*)' ' end if if (keyword.eq.'ROTPOSN') then * write(*,995)line(1:80) rotang = get_real (line, i1, i2, narg, 0.) write(*,*)rotang,' ROTPOSN = PA -90 deg' write(*,*)' ' end if if (keyword.eq.'PREPIX') then prepix = get_real (line, i1, i2, narg, 0.) * write(*,995)line(1:80) write(*,*)prepix,' Prepix pixels in CCD readout ' write(*,*)' ' end if if (keyword.eq.'PRELINE') then preline = get_real (line, i1, i2, narg, 0.) * write(*,995)line(1:80) write(*,*)preline,' Preline pixels ' write(*,*)' ' end if go to 100 190 continue close (image_lu) go to 299 C....File is NOT a FITS file 200 continue C C ************************************************ C C The section to read FITS keywords from Figaro files HAS NOT C been updated in version 3.0, and is kept for historical interest C only C write(*,*)' In section to read Figaro files ' write(*,*)' ' C C This section reads in the header parameters for a Figaro file. C call dta_asfnam('zimage',vstring(1:len),'old',0,image,istatus) if (istatus.ne.0) then call fig_dtaerr(istatus,' Error opening input file ' : //vstring(1:len)) write(*,*)' File is neither Fits nor Figaro' stop end if call dta_rdvarf('zimage.fits.rotposn',1,rotang,istatus) if (istatus.ne.0) : write(*,*)' Problem with *.fits.rotposn ',istatus call dta_rdvarf('zimage.fits.prepix',1,prepix,istatus) if (istatus.ne.0) : write(*,*)' Problem with *.fits.prepix ',istatus call dta_rdvarf('zimage.fits.preline',1,preline,istatus) if (istatus.ne.0) : write(*,*)' Problem with *.fits.preline ',istatus call dta_rdvarf('zimage.fits.equinox',1,epoch0,istatus) if (istatus.ne.0) : write(*,*)' Problem with *.fits.equinox ',istatus call dta_rdvard('zimage.fits.el',1,elevation,istatus) if (istatus.ne.0) : write(*,*)' Problem with *.fits.el ',istatus C C Now the date of the image, which is the equinox of the coordinates C call dta_rdvarc('zimage.fits.date',64,string,istatus) if (istatus.ne.0) : write(*,*)' Problem reading *.fits.date ',istatus do j=1,64 if (string(j:j) .eq. '/') string(j:j)=' ' end do nstart=1 call sla_intin(string,nstart,iday,istatus) if (istatus.ne.0) write(*,*)' Error reading day ' call sla_intin(string,nstart,imonth,istatus) if (istatus.ne.0) write(*,*)' Error reading month ' call sla_intin(string,nstart,iyear,istatus) if (istatus.ne.0) write(*,*)' Error reading year ' date=iyear+1900.0+float(imonth)/12.+float(iday)/365. if (iyear.lt.90) date=date+100. C C This seems to have changed with year 2000, and also with C LRIS red, so instead I type it in. C write(*,*)' input date of observation, year, month, date ' write(*,*)' example - 2001 5 18 ' read(*,*) iyear,imonth,iday date=iyear+float(imonth)/12.+float(iday)/365. write(*,*)' adopted date is ',date equinox=date C C Decode the declination C call dta_rdvarc('zimage.fits.dec',11,string,istatus) if (istatus.ne.0) : write(*,*)' Problem with *.fits.dec ',istatus do j=1,11 if (string(j:j) .eq. ':') string(j:j)=' ' end do do j=12,80 string(j:j)=' ' end do do j=4,80 string1(j:j)=string(j-3:j-3) end do string1(1:1)=' ' string1(2:2)=' ' string1(3:3)=' ' nstart=3 call sla_dafin(string1,nstart,dec0,istatus) if (istatus.ne.0) then write(*,*)' Error decoding *.fits.dec ',istatus end if C C Now decode the RA C nstart=1 call dta_rdvarc('zimage.fits.ra',nstart,string,istatus) if (istatus.ne.0) : write(*,*)' Problem with *.fits.ra ',istatus do j=1,11 if (string(j:j) .eq. ':') string(j:j)=' ' end do do j=12,80 string(j:j)=' ' end do do j=2,80 string1(j:j)=string(j-1:j-1) end do string1(1:1)=' ' nstart=2 call sla_dafin(string1,nstart,ra0,istatus) if (istatus.ne. 0) then write(*,*)' Error decoding *.fits.ra ',istatus end if ra0=ra0*15.0d0 call dta_rdvarc('zimage.fits.ha',11,string,istatus) if (istatus.ne.0) : write(*,*)' Problem with *.fits.ha ',istatus do j=1,11 if (string(j:j) .eq. ':') string(j:j)=' ' end do do j=12,80 string(j:j)=' ' end do do j=5,80 string1(j:j)=string(j-4:j-4) end do string1(1:4)=' ' nstart=4 call sla_dafin(string1,nstart,hour_angle,istatus) if (istatus.ne.0) then write(*,*)' Error decoding *.fits.ha ',istatus end if hour_angle=hour_angle*15.0d0 call dta_fclose('zimage',istatus) if (istatus.ne.0) then call fig_dtaerr(istatus,' Error closing image file ') end if C C End of Figaro file header section C C C ************************************************ C 299 continue write(*,*)' centerfield ra in radians (ra0) ',ra0 write(*,*)' centerfield dec in radians (dec0) ', dec0 write(*,*)' hour angle in radians ',hour_angle write(*,*)' ' C....Continue to read the input file 300 continue read (infile, lnfmt, iostat=in_stat) line if (in_stat)400,301,912 301 continue C....Parse line and determine type of input i1 = 1 i2 = 80 narg = d_parse (line, i1, i2, vstring, len) C........check for null line if (narg.eq.0) go to 300 C........check for keyword line if (vstring(1:1).eq. '\\') then keyword(1:6) = vstring(2:7) section = 0 if (keyword(1:6).eq.'REFERE') section = 2 if (keyword(1:6).eq.'OBJECT') section = 3 if (keyword(1:3).eq.'ADC') section = 4 if (keyword(1:3).eq.'END') go to 400 if (section.eq.0) go to 914 go to 300 end if C........must be data by default; reset parse; make sure data is expected i1 = 1 if (section.lt.2 .or. section.gt.4) go to 915 C........is it reference data? if (section.eq.2) then do j=1,80 if (line(j:j).eq. ':')line(j:j)=' ' end do do j=2,80 string(j:j)=line(j-1:j-1) end do string(1:1)=' ' nstart=2 call sla_flotin(string,nstart,refx(nref+1),istatus) if (istatus.ne.0) write(*,*)' Decoding problem refx ' call sla_flotin(string,nstart,refy(nref+1),istatus) if (istatus.ne.0) write(*,*)' Decoding problem refy' call sla_dafin(string,nstart,refra(nref+1),istatus) if (istatus.ne.0) write(*,*)' Decoding problem refra' call sla_dafin(string,nstart,refdec(nref+1),istatus) if (istatus.ne.0) write(*,*)' Decoding problem refdec' call sla_flotin(string,nstart,refepoch(nref+1),istatus) if (istatus.ne.0) write(*,*)' Decoding problem refepoch ' refra(nref+1)=15.0d0*refra(nref+1) if (istatus.ne.0) then write(*,*)' Problem with reference line -- skipped' go to 300 end if nref = nref + 1 if (nref. gt. max_ref) then write(*,*)' Too many reference objects, maximum is ',max_ref nref=max_ref end if endif C........or selected object? if (section.eq.3) then x(nobj+1) = get_real (line, i1, i2, narg, 0.) y(nobj+1) = get_real (line, i1, i2, narg, 0.) if (narg.eq.2) then nobj = nobj + 1 if (nobj. gt. max_obj) then write(*,*)' Too many objects, maximum is 100' nobj=max_obj end if else write (*,*)' object line format error -- skipped' endif mag(nobj) = get_real (line, i1, i2, narg, def_mag) priority(nobj) = get_int (line, i1, i2, narg, def_pri) call get_field (line, i1, i2, narg, id(nobj), 10) endif if (section.eq.4) then nstart=1 call sla_intin(line,nstart,iadc,istatus) if (istatus.ne.0) write(*,*)' Error reading ADC ' write(*,*)iadc, ' ADC status read from input file' write(*,*)' (0=ADC not in use, 1=ADC is in use)' write(*,*)' ' end if go to 300 400 continue close (infile) C C C Section for which red side Dewar was used and what scale to use in C subroutine transform. C C fit_flag=1 original Dewar, fit_flag=2 new Dewar (1996.5) C fit_flag=3 2007 fit, no ADC C fit_flag=4 2007 fit, yes ADC C C iblue=0 LRIS-R iblue=1 LRIS-B C iadc=0 No ADC (true for prior to 2007) C iadc=1 uses ADC (only possible after 2007.0) C date_install refers to new LRIS-R camera (installed in 1996) C blue_install refers to when LRIS-B mosaic CCD detector was installed C date_install = 1996.51 ! July 4, 1996 if (iblue.eq.1) then namp=4 ! 2 per CCD in the mosaic if (iadc.eq.1) then fit_flag = 10 else fit_flag=11 end if else C C Section for choice of fit for LRIS-R C if (date.gt. date_install) then fit_flag=2 else fit_flag=1 end if if (date.gt.2006.0) fit_flag=3 if (iadc.eq.1) fit_flag=4 end if write(*,*)' ' write(*,*)' ASSUMPTIONS - LRIS RED ' write(*,*)' Images taken after July 1, 1996' write(*,*)' use the new LRIS-R Dewar' write(*,*)' and images taken before then use the original' write(*,*)' LRIS-R Dewar. The scale is slightly different.' write(*,*)' ' write(*,*)' Images taken after 2006.0 use the fits derived' write(*,*)' in early 2007. Fits are available with' write(*,*)' and also without the ADC' write(*,*)' ' write(*,*)' ASSUMPTIONS - LRIS BLUE ' write(*,*)' Only fits for the mosaic CCD detector are' write(*,*)' available. Images taken with the original' write(*,*)' LRIS-B detector (a 2048x2048 CCD) cannot' write(*,*)' be processed.' write(*,*)' Fits are available with and without the ADC' write(*,*)' ' if (fit_flag.eq.1) : write(*,*)' This frame used the original LRIS-R Dewar.' if (fit_flag.eq.2) : write(*,*)' This frame uses the new LRIS-R Dewar. ' if (fit_flag.eq.3) : write(*,*)' This frame uses the fits of March 2007, no ADC' if (fit_flag.eq.4) : write(*,*)' This frame uses the fits of March 2007, ADC YES' write(*,*) if (fit_flag.eq.2 .and. date.gt. 1998.5) then write(*,*)'WARNING WARNING WARNING' write(*,*)' There was no monitoring of LRIS plate' write(*,*)' scales from late 1996 to early 2007.' write(*,*)' ' end if if (fit_flag.eq.10) write(*,*) : ' Using LRIS-B, ADC in, distortion fits of spring 2007' if (fit_flag.eq.11) write(*,*) : ' Using LRIS-B, ADC out, distortion fits of spring 2007' write(*,*) if (iadc.eq.1) write(*,*)' ADC status: in use ' if (iadc.eq.0) write(*,*)' ADC status: not in use ' write(*,*) C C End of this added section for which distortion fit to use C C C Verify input data C write(*,*)' There are ',nref,' reference stars.' write(*,*)' Their x, y, ra, dec, and epoch values are: ' write(*,*)' (RA and DEC are in radians now.) ' write(*,*)' ' write (*,105)(refx(i), refy(i), refra(i), refdec(i), : refepoch(i), i=1, nref) 105 format('ref: ',2f7.1,2f9.4,f7.1) write(*,*)' ' write(*,*)' There are ',nobj,' object stars.' write(*,*)' Their x, y, ra, mag, priority, and id values are: ' write(*,*)' ' write (*,106) (x(i), y(i), mag(i), priority(i), id(i), i=1, nobj) 106 format('obj: ',2f7.1,f6.2,I7,X,A10) write(*,*)' The equinox of the object data is ',equinox C............................................................................. C First step - subtract the pre-image pixels (2*PREPIX) from each x C and the preline pixels (PRELINE) from each y. C do i=1,nref refx(i)=refx(i)-float(namp)*prepix refy(i)=refy(i)-preline end do do i=1,nobj x(i)=x(i)- float(namp)*prepix y(i)=y(i)-preline end do value=float(namp)*prepix write(*,*)' ' write(*,*)' Input data corrected for ',value,' pre-image pixels' write(*,*)' Input data corrected for ',preline,' pre-image lines' write(*,*)' ' C C Next apply transformation from (x,y) in pixels to xarc,yarc in C arc-seconds. Details are given in the subroutine documentation C The (xarc,yarc) coordinate system has positive xarc to the right, C positive yarc is up, i.e. a standard right handed coordinate system. C open(unit=14, file = 'arcsec.list') if (iblue.eq.0) then do i=1,nref call transform_red(refx(i),refy(i),xarc_ref(i),yarc_ref(i), : fit_flag) write(14,*)i,refx(i),refy(i),xarc_ref(i),yarc_ref(i), : ' Transform R REF' end do do i=1,nobj call transform_red(x(i),y(i),xarc(i),yarc(i), : fit_flag) end do write(*,*)' Coordinates transformed from pixels to arc-sec.' write(*,*)' LRIS-Red ' write(*,*)' ' else do i=1,nref call transform_blue(refx(i),refy(i),xarc_ref(i),yarc_ref(i), : fit_flag) write(14,*)i,refx(i),refy(i),xarc_ref(i),yarc_ref(i), : ' Transform B REF' end do write(14,*)' ' do i=1,nobj call transform_blue(x(i),y(i),xarc(i),yarc(i), : fit_flag) write(14,*)i,x(i),y(i),xarc(i),yarc(i),' Transform B OBJ' end do write(*,*)' Coordinates transformed from pixels to arc-sec.' write(*,*)' LRIS-Blue ' write(14,*)' LRIS-Blue used here' write(*,*)' ' end if C C Now we must rotate so that (xarc,yarc) becomes (north_arc,east_arc) C Note that the position angle in conventional astronomical usage C is the keyword returned by the Keck drive and control system, C ROTPOSN (given in degrees) + 90.0 degrees. Note that C East is -x. C angle=rotang+90.0 angle=-angle C C Enter here for 2nd try with new PA from reference stars C 320 continue C C Correction of bug, April 18, 1995 C do i=1,max_ref value=cosd(angle)*xarc_ref(i) +sind(angle)*yarc_ref(i) east_arc_ref(i)=-value north_arc_ref(i)= -sind(angle)*xarc_ref(i) : +cosd(angle)*yarc_ref(i) end do do i=1,max_obj value=cosd(angle)*xarc(i) +sind(angle)*yarc(i) east_arc_obj(i)=-value north_arc_obj(i)= -sind(angle)*xarc(i)+cosd(angle)*yarc(i) end do write(*,*)' Rotation correction applied, position angle = ', : angle write(*,*)' ' write(14,*)' Reference stars: delta(RA), delta(Dec) (arc-sec)' write(14,*)' with respect to the center of the field.' write(14,*)' No atmospheric refraction correction applied.' write(14,*)' ' do i=1,nref write(14,550)i,east_arc_ref(i),north_arc_ref(i) 550 format(1x,i3,3x,f7.2,2x,f7.2) end do write(14,*)' ' write(14,*)' Objects now: ' write(14,*)' ' do i=1,nobj write(14,551)i,east_arc_obj(i),north_arc_obj(i),id(i) 551 format(1x,i3,2x,f7.2,2x,f7.2,5x,a10) end do write(*,*)' Arc sec N and E of center offsets in file' write(*,*)' arcsec.list' C C Now find the correction for atmospheric refraction for each point. C This is the correction to turn apparent RA, dec into true RA, dec. C C First compute the apparent RA and DEC for each point using the C centerfield coordinates, ra0, dec0 (in radians), and adding the C the offsets. Don't forget cos(dec) for RA. C Then call subroutine refraction. rpas = radians/arc-sec C rpas=pi/(180.0*3600.0) * dvalue= dcos(dec0) do i=1,nref adec_ref(i)=dec0 +north_arc_ref(i)*rpas dvalue= dcos(adec_ref(i)) ara_ref(i)=ra0 + east_arc_ref(i)*rpas/dvalue dvalue1=hour_angle - east_arc_ref(i)*rpas/dvalue call refraction(ara_ref(i),adec_ref(i),true_ra_ref(i), : true_dec_ref(i),elevation,dvalue1) end do do i=1,nobj adec_obj(i)=dec0 + north_arc_obj(i)*rpas ara_obj(i)=ra0 + east_arc_obj(i)*rpas/dvalue dvalue= dcos(adec_obj(i)) dvalue1=hour_angle - east_arc_obj(i)*rpas/dvalue call refraction(ara_obj(i),adec_obj(i),true_ra_obj(i), : true_dec_obj(i),elevation,dvalue1) end do C C refraction for the centerfield coordinates C call refraction(ra0,dec0,true_ra_center,true_dec_center, : elevation, hour_angle) write(*,*)' Refraction corrections applied.' write(*,*)' ' C C Just as a check for the astronomer, calculate the maximum and C minimum refraction in RA and DEC in arc-sec with respect to the C center of the field. C do i=1,nref delta_ra(i)=(true_ra_ref(i) -ara_ref(i) - (true_ra_center : -ra0))/rpas delta_dec(i)=(true_dec_ref(i) -adec_ref(i) -(true_dec_center : -dec0))/rpas end do dra_max=-100.0d0 dra_min=100.0d0 ddec_max=-100.0d0 ddec_min=100.0d0 do i=1,nref if (delta_ra(i) .gt. dra_max) dra_max=delta_ra(i) if (delta_ra(i ).lt. dra_min) dra_min=delta_ra(i) if (delta_dec(i) .gt. ddec_max) ddec_max=delta_dec(i) if (delta_dec(i) .lt. ddec_min) ddec_min=delta_dec(i) end do write(*,*)' For reference stars: ' write(*,*)' Maximum and minimum refraction correction in arc-sec' write(*,*)' compared to the center field for RA ' write(*,*)dra_max,dra_min write(*,*)' ' write(*,*)' Maximum and minimum refraction correction in arc-sec' write(*,*)' compared to the center field for DEC ' write(*,*)ddec_max, ddec_min write(*,*)' ' write(*,*)' ' C C Now for object list C do i=1,nobj delta_ra(i)=(true_ra_obj(i)-ara_obj(i) - (true_ra_center-ra0)) : /rpas delta_dec(i)=(true_dec_obj(i)-adec_obj(i)- : (true_dec_center-dec0))/rpas end do dra_max=-100.0d0 dra_min=100.0d0 ddec_max=-100.0d0 ddec_min=100.0d0 do i=1,nobj if (delta_ra(i) .gt. dra_max) dra_max=delta_ra(i) if (delta_ra(i ).lt. dra_min) dra_min=delta_ra(i) if (delta_dec(i) .gt. ddec_max) ddec_max=delta_dec(i) if (delta_dec(i) .lt. ddec_min) ddec_min=delta_dec(i) end do write(*,*)' Maximum and minimum refraction correction in arc-sec' write(*,*)' compared to the center field for RA for objects ' write(*,*)dra_max, dra_min write(*,*)' ' write(*,*)' Maximum and minimum refraction correction in arc-sec' write(*,*)' compared to the center field for DEC for objects ' write(*,*)ddec_max, ddec_min write(*,*)' ' C C Calculate final coordinates. If nref = 0, then use the C centerfield values - they already include a refraction correction, C put in by the Keck DCS, C I believe, so that has to be re-zeroed to the centerfied value. C C If nref > 0, use average offset. Print out each offset so C user can check. Also print out mean and sigma of the offsets. C User is allowed to reject one and only one reference star. C C if (nref .eq. 0) then do i=1,nobj true_ra_obj(i)= ra0 + (true_ra_obj(i)-true_ra_center) true_dec_obj(i) = dec0 + (true_dec_obj(i)-true_dec_center) end do write(*,*)' Object coordinates set to centerfield position' write(*,*)' ' go to 600 end if C C First try using all the reference stars C nref_bad=0 601 continue j=1 do i=1,nref if (ref_ok(i).eq.0) then delta_ra(j)=true_ra_ref(i)-refra(i) delta_dec(j)=true_dec_ref(i)-refdec(i) j=j+1 end if end do j=j-1 write(*,*)' Results of reference star offsets in arc-sec' write(14,*)' Results of reference star offsets in arc-sec' j=1 do i=1,nref if (ref_ok(i).eq.0) then dvalue=delta_ra(j)/rpas dvalue1=delta_dec(j)/rpas write(*,140)i, dvalue,dvalue1 write(14,140)i, dvalue,dvalue1 140 format(1x,i3,' Delta RA: ',e16.8,1x,' Delta DEC: ',e16.8) j=j+1 else write(*,*)' Reference star ',i,' was eliminated.' end if end do j=j-1 write(*,*)' ' if (j.gt.1) then call sigav(delta_ra,j,ra_average,sigma_ra) call sigav(delta_dec,j,dec_average,sigma_dec) write(*,*)' There are ',j,' good reference stars.' write(*,*)' ' dvalue=ra_average/rpas dvalue1=sigma_ra/rpas write(*,141)dvalue,dvalue1 141 format(' Mean Delta RA, and sigma: ',2e16.8) dvalue=dec_average/rpas dvalue1=sigma_dec/rpas write(*,142)dvalue, dvalue1 142 format(' Mean Delta DEC, and sigma: ',2e16.8) write(*,*)' ' end if write(*,*) ' You now have 3 choices: ' write(*,*)' Choice 1: Use this solution. Type 1.' write(*,*) : ' Choice 2: Eliminate a reference star and re-solve. Type 2. ' write(*,*)' Choice 3: Use the centerfield only. Type 3 .' 602 read(*,*)ichoice if (ichoice .ne. 1 .and. ichoice.ne.2 .and. ichoice.ne.3) then write(*,*)' You have not chosen either choice 1, 2, or 3.' write(*,*)' Please indicate your choice again. ' go to 602 end if if (ichoice .eq. 3) then do i=1,nobj true_ra_obj(i)= ra0 + (true_ra_obj(i)-true_ra_center) true_dec_obj(i) = dec0 + (true_dec_obj(i)-true_dec_center) end do write(*,*)' Object coordinates set to centerfield position' write(*,*)' ' go to 600 end if if (ichoice .eq. 1) then do i=1,nobj true_ra_obj(i)=true_ra_obj(i) - ra_average true_dec_obj(i)=true_dec_obj(i) - dec_average end do write(*,*)' Object coordinates set to reference stars' write(*,*)' ' go to 600 end if if (ichoice.eq.2) then if (nref_bad.gt.0) then write(*,*)' You have already eliminated 1 reference star.' write(*,*)' You cannot eliminate more than 1.' write(*,*)' Please choose again, options 1 or 3 only. ' go to 602 end if if (nref .lt. 4) then write(*,*)' You have less than 4 reference stars.' write(*,*)' You are not allowed to eliminate any. ' write(*,*)' Please choose again, options 1 or 3 only.' go to 602 end if write(*,*)' What is the number of the reference star you' write(*,*)' wish to eliminate ? ' read(*,*)nref_bad ref_ok(nref_bad)=100 go to 601 end if 600 continue C C If there are 10 or more reference stars, we compute the PA from C the set of reference stars. We do this using C /scr/jlc/keckastrom/matrix.f followed by test_rotate.f C C we use xarc_ref(i),yarc_ref(i) and xxx(i),yyy(i) as the input C to the matrix solver. C if (nref.ge.10 ) then ! big loop for solving PA from reference stars do i=1,nref dvalue = dcos(refdec(i)) xtrue(i) = xarc_ref(i) ytrue(i) = yarc_ref(i) xxx(i)= dvalue*(refra(i) - (ra0 - ra_average))/rpas yyy(i)= (refdec(i) - (dec0 - dec_average))/rpas end do C C solve for centerfield location and position angle from objects C on LRIS image C call solve_matrix(nref,xtrue,ytrue,xxx,yyy, : deltax,deltay,theta_image) write(*,*)' ' write(*,192)deltax,deltay,theta_image 192 format('delta X,Y of centerfield ',2(2x,f12.3),5x, : ' Mean angle ',f12.3) C C Compare theta from header with theta from image C theta_image = theta_image -180.0 !!!CHECK for ALL frames!!! write(*,193)angle,theta_image 193 format(2(2x,f12.3),' PA from header, PA from image ') if (abs(angle-theta_image).gt. pa_tol) then write(*,*)' ' write(*,*)' Difference (deg) bigger than ',pa_tol write(*,*)' Do you want to fix this ?' write(*,*)' 0 = no, 1 = yes ' read(*,*)icheck if (icheck.eq.1) then angle = theta_image write(*,*)' Position angle changed to ',angle go to 320 end if end if C C ++++ end of PA check section, need to print out, ask if implement...NOT DONE YET C end if close(unit=14) write(*,*)' ' C C Now convert final RA,DEC into strings, and write output file C Suitable for use with Autoslit. C C SLALIB routines used here: C sla_dr2tf, radians to hours/minutes/seconds C sla_dr2af, radians to degrees/ arc-minutes / arc-seconds C do i=1,nobj nstart=4 call sla_dr2tf(nstart,true_ra_obj(i),cvalue_ra(i),int_array) do j=1,4 iarray_ra(j,i)=int_array(j) end do call sla_dr2af(nstart,true_dec_obj(i),cvalue_dec(i), : int_array) do j=1,4 iarray_dec(J,i)=int_array(j) end do end do C C Write the output file for Autoslit. First step is to fix up C writing the sign of the DEC. The format I4.4 puts in leading C zeros as appropriate for the decimal fraction the RA and DEC C seconds (of arc and of time). The first line is for the C center of the field. We use the average RA and DEC for the C set of objects. C dsum1=0.0d0 dsum2=0.0d0 do i=1,nobj dsum1=dsum1+true_ra_obj(i) dsum2=dsum2+true_dec_obj(i) end do center_obj_ra=dsum1/float(nobj) center_obj_dec=dsum2/float(nobj) call sla_dr2tf(nstart,center_obj_ra,cvalue_center_ra,icenter_ra) call sla_dr2af(nstart,center_obj_dec,cvalue_center_dec, : icenter_dec) C C Now we are ready to write the output file. C open(unit=12,file='coordinates.list') if (cvalue_center_ra .eq. '-') : write(*,*)' Error in ra calculation for centerfield. ' C C Original 5 lines above here replaced by a section with the proper treatment C of the sign of the declination. (July 1997). Same thing C several lines down to output the declination of the object list. C dchar='+' if (cvalue_center_dec .eq. '-') dchar='-' encode(2,149,string2)icenter_dec(1) 149 format(i2) string3(1:1)=dchar(1:1) string3(2:3)=string2(1:2) if (icenter_dec(1) .lt. 10.0) string3(2:2)='0' write(12,143)(icenter_ra(i),i=1,4),string3(1:3), : (icenter_dec(i),i=2,4),epoch0, : equinox 143 format(' CENTER 9999 0.00 ',i3,1x,i3,1x,i3,'.',i4.4,3x,a3,1x, : i3,1x,i3,'.',i4.4,2x,f7.1,2x,f7.1,' 0.0 0.0 ') do i=1,nobj if (cvalue_ra(i) .eq. '-') : write(*,*)' Error in ra calculation.' dchar='+' if (cvalue_dec(i) .eq. '-') dchar='-' encode(2,149,string2)iarray_dec(1,i) string3(1:1)=dchar(1:1) string3(2:3)=string2(1:2) if (iarray_dec(1,i) .lt. 10.0) string3(2:2)='0' write(12,144)id(i),priority(i),mag(i),(iarray_ra(j,i), : j=1,4),string3(1:3), : (iarray_dec(j,i),j=2,4), refepoch(1),equinox 144 format(a10,2x,i5,2x,f7.2,3x,i3,1x,i3,1x,i3,'.',i4.4, : 3x,a3,1x,i3,1x,i3,'.',i4.4,2x,f7.1,2x,f7.1, : ' 0.0 0.0 ') end do close(unit=12) write(*,*)' Final output is file coordinates.list' write(*,*)' ' stop 910 write (*,*)'Error on opening input file (does not exist?)' stop 911 write (*,*)'Error on opening image file' stop 912 write (*,*)'Error reading input file' stop 913 write (*,700)in_stat 700 format(' Error reading image header',I6) stop 914 write (*,701)vstring(1:len) 701 format(A20,' : Unknown keyword!') stop 915 write (*,702)line(1:80) 702 format(A80,' Incorrect input format! (missing keyword?)') stop 916 write (*,*)' Image file does not exist! ' stop end C NOTES: c ideally, should allow input line-length to change easily (set at 80 char) c [FYI] Only the first six characters of the keyword are checked (by design) C+ subroutine sigav (array,index,average,sigma) C C Computes mean and standard deviation for an array with index C entries of double precision numbers. array(index) C real*8 array(index),sum,average,sigma,value sum=0.0 do i=1,index sum=sum+array(i) end do average=sum/float(index) sum=0.0 do i=1,index sum=sum+(array(i)-average)*(array(i)-average) end do value = sum/(float(index)-1.0) sigma=dsqrt(value) return end C+ subroutine transform_red(xin,yin,xout,yout,fit_flag) C C for LRIS-R C Applies the transformation from (x,y) in pixels on the CCD with C the origin of the coordinates at the lower left corner of the C CCD, so the prescan and preline pixels are not counted. C C C The output is x,y in arc-sec with the origin at the center of the C 2048x2048 CCD. C C Fits 1 and 2 are from "The Spatial Distortion in the LRIS" C by J.Cohen, Sep. 12, 1994, page 31. C C The output coordinates (xout,yout) have a standard right hand C coordinate system with x increasing toward the right and y C increasing upward. C C Modified Sep. 1996 C Added parameter "fit_flag" (an integer) C This parameter takes care of the difference in scale between C the original LRIS dewar and the Dewar installed in July 1996 C with the HIRES CCD and the meniscous field flattener designed C by Jim McCarthy. C implicit none real*8 a(10),b(10),c(10),d(10),e(10),f(10),g(10),h(10) real*8 gno(10), hno(10), xsq,ysq real*4 xin,yin,xout,yout,scale integer fit_flag,fit_check C C These are the fits for the original LRIS dewar. (fit_flag=1) C data a /-0.21727366E+03, 0.21088950E+00, -0.36604795E-02, : 0.29435336E-05, 0.18486891E-05, 0.25794732E-05, : -0.11709022E-08, 0.46571952E-10, : -0.13710488E-08, 0.11776678E-10 / data b / 0.21765796E+03, 0.17362377E-02, -0.20996992E+00, : -0.16847447E-05, -0.37531701E-05, -0.15077788E-05, : 0.12859979E-09, 0.12886424E-08, -0.16780111E-10, : 0.12571795E-08 / C C These are the fits for the new LRIS dewar with the meniscous field C flattener installed July 1996. The difference in scale is about C a factor of 1.0105. The coefficients are those of the M71 frame C taken in Sep 1996. (fit_flag=2) C data c / 0.22218600E+03, -0.20837592E+00, 0.40227488E-02, : -0.29374097E-05, -0.20568340E-05, -0.23432073E-05, : 0.10152123E-08, 0.85416046E-10, 0.10912104E-08, : 0.10591175E-09 / data d / -0.21475296E+03, -0.12604107E-02, 0.20844918E+00, : 0.10091565E-05, 0.27702410E-05, 0.15872904E-05, : 0.45784514E-10, -0.10796130E-08, -0.12455420E-09, : -0.84843272E-09 / C C Fits to LRIS-R, YES ADC, done from frame lred0040 of NGC 2158 C taken in March 2007 C data g / 0.21539475E+03, -0.20950212E+00, 0.30710551E-02, : -0.20525413E-05, : -0.14337209E-05, -0.21997301E-05, 0.77741210E-09, : -0.17289732E-10, 0.11033589E-08, -0.57206933E-10 / data h / -0.21610848E+03, -0.11223419E-02, 0.20869838E+00, : 0.11890912E-05, : 0.33875056E-05, 0.10158390E-05, -0.57186568E-10, : -0.10449538E-08, 0.48144598E-10, -0.11278520E-08 / C C Fits to LRIS-R, NO ADC, done from frame lred0070 of NGC 2158 C taken in March 2007 C data gno / 0.21483418E+03, -0.20825571E+00, 0.32599909E-02, : -0.30204635E-05, : -0.15392512E-05, -0.23929830E-05, 0.10442439E-08 , : 0.11155732E-09, 0.10654012E-08 , -0.14505548E-10 / data hno / -0.21568238E+03, -0.11120534E-02, 0.20829835E+00, : 0.11878170E-05, : 0.32459433E-05, 0.11281612E-05, -0.68787117E-10, : -0.10104519E-08, -0.70793360E-11, -0.10475806E-08 / C C C fit_flag=1 original Dewar, C fit_flag=2 new Dewar, 1996 upgrade C fit_flag=3 2007, LRIS-R, no ADC C fit_flag=4 2007, LRIS-R, yes ADC C xsq=xin*xin ysq=yin*yin fit_check=0 if (fit_flag.eq.1) then fit_check=1 xout=a(1) +a(2)*xin + a(3)*yin +a(4)*xsq + a(5)*ysq : +a(6)*xin*yin +a(7)*xin*xsq +a(8)*xsq*yin : +a(9)*xin*ysq + a(10)*yin*ysq yout=b(1)+b(2)*xin +b(3)*yin +b(4)*xsq +b(5)*ysq : +b(6)*xin*yin +b(7)*xin*xsq +b(8)*xsq*yin : +b(9)*xin*ysq + b(10)*yin*ysq C C Now this gives +xout is increasing to the right, so OK C and +yout is increasing downward. We must flip yout. C yout = -yout end if if (fit_flag.eq.2) then fit_check=1 xout=c(1) +c(2)*xin + c(3)*yin +c(4)*xsq + c(5)*ysq : +c(6)*xin*yin +c(7)*xin*xsq +c(8)*xsq*yin : +c(9)*xin*ysq + c(10)*yin*ysq yout=d(1)+d(2)*xin +d(3)*yin +d(4)*xsq +d(5)*ysq : +d(6)*xin*yin +d(7)*xin*xsq +d(8)*xsq*yin : +d(9)*xin*ysq + d(10)*yin*ysq C C with this set of coefficients x but not y must be flipped C xout=-xout end if if (fit_flag.eq.3) then fit_check=1 xout=gno(1) +gno(2)*xin + gno(3)*yin +gno(4)*xsq + : gno(5)*ysq : +gno(6)*xin*yin +gno(7)*xin*xsq +gno(8)*xsq*yin : +gno(9)*xin*ysq + gno(10)*yin*ysq yout=hno(1)+hno(2)*xin +hno(3)*yin +hno(4)*xsq : +hno(5)*ysq : +hno(6)*xin*yin +hno(7)*xin*xsq +hno(8)*xsq*yin : +hno(9)*xin*ysq + hno(10)*yin*ysq C C with this set of coefficients x but not y must be flipped C xout=-xout end if if (fit_flag.eq.4) then fit_check=1 xout=g(1) +g(2)*xin + g(3)*yin +g(4)*xsq + g(5)*ysq : +g(6)*xin*yin +g(7)*xin*xsq +g(8)*xsq*yin : +g(9)*xin*ysq + g(10)*yin*ysq yout=h(1)+h(2)*xin +h(3)*yin +h(4)*xsq +h(5)*ysq : +h(6)*xin*yin +h(7)*xin*xsq +h(8)*xsq*yin : +h(9)*xin*ysq + h(10)*yin*ysq C C with this set of coefficients x but not y must be flipped C xout=-xout end if if (fit_check.ne.1) then write(*,*)' Error - scale not determined in' write(*,*)' subroutine transform. ',fit_flag stop end if return end C+ subroutine transform_blue(xin,yin,xout,yout,fit_flag) C C for LRIS-Blue C Applies the transformation from (x,y) in pixels on the CCD with C the origin of the coordinates at the lower left corner of the C CCD, so the prescan and preline pixels are not counted. C C C The output is x,y in arc-sec with the origin at the center of the C 2048x2048 CCD. C C The fits used are from frames lblue0055 and lblue0062, ADC C in and out, from spring 2007. See files C /scr/jlc/keckastrom/adc_april2007/notes* C C The output coordinates (xout,yout) have a standard right hand C coordinate system with x increasing toward the right and y C increasing upward. C C implicit none real*8 a(10),b(10),c(10),d(10),xsq,ysq real*8 xnew, ynew real*4 xin,yin,xout,yout,scale,xstart,ystart real*4 xcon, ycon real*4 gap_x, gap_y, theta, delx, dely, aval integer fit_flag, fit_check, iend_left C C These are the fits for LTIS-B with the ADC IN C data a / -0.27483492E+03, -0.10693314E-02, 0.13192481E+00, : 0.24581167E-06, : 0.12822328E-05, 0.84859900E-06, -0.43086048E-11, : -0.19709703E-09, -0.35997436E-11, -0.15949522E-09 / data b / -0.27976069E+03, 0.13151461E+00, -0.22063626E-02, : 0.12483489E-05, 0.25379359E-06 , 0.11612311E-05, : -0.19916774E-09, -0.53179088E-12, : -0.16899998E-09, 0.14019297E-10 / C C These are the fits LRIS-B with the ADC OUT C data c / -0.27442599E+03, -0.12496566E-02, 0.13181624E+00, : 0.32606807E-06, : 0.12004263E-05, 0.87444933E-06, -0.10594879E-10, : -0.20943796E-09, 0.21221213E-11, -0.14837532E-09 / data d / -0.27870945E+03, 0.13102987E+00, -0.26848773E-02, : 0.12898216E-05, : 0.40423847E-06, 0.13085340E-05, -0.20381190E-09, : -0.12707996E-11, -0.20420376E-09, 0.65189775E-12 / C C Constants for CCD mosaic geometry C iend_left = 2048 aval = float(iend_left) xcon = 2896.0 ! offsets to get to center of right side ycon = 2125.0 theta = -0.092 ! rotation between 2 CCDs in mosaic gap_x = 99.92 ! gap between 2 CCDs in mosaic gap_y = -4.43 C C C fit_flag=10 Blue side, ADC IN, fit_flag=11 BLUESIDE ADC OUT C C C First step is to transform the coordinates to take into C account the geometry of the LRIS-B CCD mosaic. There is C rotation between the 2 CCDs (a small amount) and a gap between them. C Wenjin figured all this out. C if (xin.gt. aval) then delx = xin - xcon dely = yin - ycon xnew = delx*cosd(theta) - dely*sind(theta) ynew = delx*sind(theta) + dely*cosd(theta) xstart = xnew + xcon + gap_x ystart = ynew + ycon + gap_y * * xds9(i) = xout + float(nprepix) ! nprepix = 4xprepix = 4x51 = 204 * else xstart = xin ystart = yin * xds9(i) = xout + float(nprepix) end if C C Now we apply the transformations C xsq=xstart*xstart ysq=ystart*ystart fit_check=0 if (fit_flag.eq.10) then fit_check=1 xout=a(1) +a(2)*xstart + a(3)*ystart +a(4)*xsq : + a(5)*ysq : +a(6)*xstart*ystart +a(7)*xstart*xsq + : a(8)*xsq*ystart : +a(9)*xstart*ysq + a(10)*ystart*ysq yout=b(1)+b(2)*xstart +b(3)*ystart + : b(4)*xsq +b(5)*ysq : +b(6)*xstart*ystart +b(7)*xstart*xsq + : b(8)*xsq*ystart : +b(9)*xstart*ysq + b(10)*ystart*ysq C C We must flip xout ? C xout = -xout end if if (fit_flag.eq.11) then fit_check=1 xout=c(1) +c(2)*xstart + c(3)*ystart +c(4)*xsq : + c(5)*ysq : +c(6)*xstart*ystart +c(7)*xstart*xsq + : c(8)*xsq*ystart : +c(9)*xstart*ysq + c(10)*ystart*ysq yout=d(1)+d(2)*xstart +d(3)*ystart + : d(4)*xsq +d(5)*ysq : +d(6)*xstart*ystart +d(7)*xstart*xsq + : d(8)*xsq*ystart : +d(9)*xstart*ysq + d(10)*ystart*ysq C C with this set of coefficients x but not y must be flipped ? C xout=-xout end if if (fit_check.ne.1) then write(*,*)' Error - scale not determined in' write(*,*)' subroutine transform.' stop end if return end C+ C+ subroutine solve_matrix(ncount,sky_ra,sky_dec,xccd,yccd, : deltax,deltay,theta) C C This code is from the program /scr/jlc/keckastrom/matrix.f C sky_ra,sky_dec(npt) RA,Dec offsets from centerfield in arcsec C xccd,yccd(npt) are X,Y(CCD) run through the distortion correction C 3rd order polynomial and turned into X,Y offsets in arcsec C output is deltax,deltay and theta C C C The relevant equations are: C C xtrue = deltax + x*cos(theta) + y*sin(theta) C ytrue = deltay - x*sin(theta) + y*cos(theta). C C IFLAG is an error flag. If it is not zero, an error has occured C during the solution. C C ****************************** C Feb 2008 - Version 3.1 has 3 additional lines towards end C of this subroutine to fix up problem with PA when C PA is in certain quadrants. See directory C /home/gringo/jlc/lris/coordinates/april2007/complaint C and also file C /home/gringo/jlc/lris/coordinates/april2007/test_images/testing_feb2008 C C ******************************* C C parameter (nmax=1000) real*8 pi, two_pi real*8 sky_ra(nmax),sky_dec(nmax) real*8 ax(10,10),dvalue(10),dummy(10),ca(10),cb(10),xv(10) real*8 yv(10), xccd(nmax),yccd(nmax),xin(nmax),yin(nmax) real*8 xout(nmax),yout(nmax) real*8 dx(nmax), dy(nmax),xaverage,yaverage,sigma C C Set arrays C ncoef = 3 * write(*,*)ncount,' in SOLVE MATRXI Ncount ' do i=1,ncount xin(i)=sky_ra(i) yin(i)= sky_dec(i) xout(i) = xccd(i) yout(i) = yccd(i) end do * write(*,*)' xin, yin, xout, yout ' * do i=1,ncount * write(*,*)xin(i),yin(i),xout(i),yout(i) * end do C C Initialize matrices C do i=1,10 xv(i)=0.0d0 yv(i)=0.0d0 do j=1,10 ax(i,j)=0.0d0 end do end do ncoef = 3 ! linear fit do n=1,ncount dvalue(1)=1.0d0 dvalue(2)=xin(n) dvalue(3)=yin(n) dvalue(4)=xin(n)*xin(n) dvalue(5)=yin(n)*yin(n) dvalue(6)=xin(n)*yin(n) dvalue(7)=dvalue(4)*dvalue(2) dvalue(8)=dvalue(4)*yin(n) dvalue(9)=xin(n)*dvalue(5) dvalue(10)=yin(n)*dvalue(5) do i=1,10 xv(i)=xv(i)+xout(n)*dvalue(i) yv(i)=yv(i)+yout(n)*dvalue(i) end do do i=1,10 ax(1,i)=ax(1,i) + dvalue(i) end do do i=1,10 ax(2,i) = ax(2,i) +dvalue(2)*dvalue(i) end do do i=1,10 ax(3,i) = ax(3,i) +dvalue(3)*dvalue(i) end do do i=1,10 ax(4,i) = ax(4,i) +dvalue(4)*dvalue(i) end do do i=1,10 ax(5,i) = ax(5,i) +dvalue(5)*dvalue(i) end do do i=1,10 ax(6,i) = ax(6,i) +dvalue(6)*dvalue(i) end do do i=1,10 ax(7,i) = ax(7,i) +dvalue(7)*dvalue(i) end do do i=1,10 ax(8,i) = ax(8,i) +dvalue(8)*dvalue(i) end do do i=1,10 ax(9,i) = ax(9,i) +dvalue(9)*dvalue(i) end do do i=1,10 ax(10,i) = ax(10,i) +dvalue(10)*dvalue(i) end do end do C C Call matrix inversion routine C call dmativ(ax,dummy,ncoef,10,sout) C C Answer C do i=1,ncoef ca(i)=0.0d0 cb(i)=0.0d0 do j=1,ncoef ca(i)=ca(i) + ax(i,j)*xv(j) cb(i)=cb(i) + ax(i,j)*yv(j) end do end do write(14,101)(ca(i),i=1,10) 101 format(2x,4(e16.8,4x)) write(14,*)' ' write(14,101)(cb(i),i=1,10) C C Calculate the residuals - line with ~/subs/sigav.f C do i=1,ncount xval=ca(1)+ca(2)*xin(i)+ca(3)*yin(i) dx(i)=xout(i)-xval yval=cb(1)+cb(2)*xin(i)+cb(3)*yin(i) dy(i)=yout(i)-yval * write(*,*)i,xin(i),yin(i),xout(i),yout(I),dx(i),dy(i) write(14,110)i,xin(i),yin(i),xout(i),yout(I),dx(i),dy(i) end do call sigav(dx,ncount,xaverage,sigma) write(*,104)ncount write(14,101)ncount 104 format(1x,' No. of points ',i4) write(*,102)xaverage,sigma write(14,102)xaverage,sigma 102 format(1x,' X deviations: average, sigma ',2x,e16.8,2x,f8.2) call sigav(dy,ncount,yaverage,sigma) write(*,103)yaverage,sigma write(14,103)yaverage,sigma 103 format(1x,' Y deviations: average sigma ',2x,e16.8,2x,f8.2) 110 format(i4,2x,2f10.3,5x,2f10.3,5x,2f10.3) C C Calculate the position angle C C C C Calculate rotation angle (code from C /scr/jlc/keckastrom/test_rotate.f) C C *********************this section of code was modified to C have the correct behavior in the four quadrants 0-90,90-180...deg C Feb 2008 C C C ca(3) = sind(angle), ca(2) = cosd(angle) C -cb(2) = sind(angle), cb(3) = cosd(angle) C C The next 2 lines were the entire original code. C a1 = atand(ca(3)/ca(2) ) a2 = atand(-cb(2)/cb(3) ) C C average of two independent values for angle C as in original code C theta = (a1 + a2)/2.0 C C Fix to get the correct quadrant, section added Feb 2008 C if (ca(3).gt. 0.0 .and. ca(2).lt. 0.0) theta = theta+180.0 if (ca(3).lt. 0.0 .and. ca(2).gt. 0.0) theta = theta+180.0 write(*,111)a1,a2,theta write(14,112)a1,a2,ca(3),ca(2),-cb(2),cb(3) ! checking quadrant write(14,111)a1,a2,theta 111 format(1x,'2 values for angle, mean ',f9.4,5x,2(f9.4,2x)) 112 format(3(2f9.3,5x),' TEST: internal angles, a1, a2 - sin,cos') deltax = ca(1) deltay = cb(1) return end C C+ C SUBROUTINE DMATIV(ARRAY,VECTOR,MORD,NORDER,DET) C+ C VERSION A 10 AUG 83 KOO,DAVID DTM VAX C C SUBROUTINE TO EMULATE DMATIV OF G. BRUZUAL C USING MATINV FROM BEVINGTON C- REAL*8 ARRAY(NORDER,NORDER) REAL*8 BRRAY(10,10),VECTOR(1),WECTOR(10) C FIRST LOAD BRRAY WITH CONTENTS OF ARRAY DO 1 I=1,MORD DO 2 J=1,MORD 2 BRRAY(J,I)=ARRAY(J,I) 1 CONTINUE C CALL MATINV(BRRAY,MORD,DET) C C NOW RELOAD OUTPUT ARRAY DO 11 I=1,MORD DO 12 J=1,MORD 12 ARRAY(J,I)=BRRAY(J,I) 11 CONTINUE C C NOW GET OUTPUT VECTOR DO 10 I=1,MORD WECTOR(I)=0. DO 20 J=1,MORD WECTOR(I)=WECTOR(I)+ARRAY(J,I)*VECTOR(J) 20 CONTINUE 10 CONTINUE DO 30 I=1,MORD 30 VECTOR(I)=WECTOR(I) RETURN END C SUBROUTINE MATINV(ARRAY,NORDER,DET) C+ C VERSION B 10 AUG 83 ADAPTED FOR VAX C VERSION A 24 MAR 80 KOO,DAVID C FROM PAGE 302 OF BEVINGTON BOOK(1969) C C PURPOSE:INVERT A SYMMETRIC MATRIX AND GET DETERMINANT C USAGE: CALL MATINV(ARRAY,NORDER,DET) C PARAMETERS: ARRAY=INPUT MATRIX WHICH IS REPLACED BY INVERSE C NORDER =DEGREE OF MATRIX(ORDER OF DETERM) C DET =DETERMINANT OF INPUT MATRIX C MAX FOR NORDER=10 C- REAL*8 ARRAY,AMAX,SAVE DIMENSION ARRAY(10,10),IK(10),JK(10) 10 DET=1. 11 DO 100 K=1,NORDER C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX AMAX=0. 21 DO 30 I=K,NORDER DO 30 J=K,NORDER 23 IF(DABS(AMAX)-DABS(ARRAY(I,J))) 24,24,30 24 AMAX=ARRAY(I,J) IK(K)=I JK(K)=J 30 CONTINUE C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K) 31 IF(AMAX) 41,32,41 32 DET=0. GO TO 140 41 I=IK(K) IF(I-K) 21,51,43 43 DO 50 J=1,NORDER SAVE=ARRAY(K,J) ARRAY(K,J)=ARRAY(I,J) 50 ARRAY(I,J)=-SAVE 51 J=JK(K) IF(J-K) 21,61,53 53 DO 60 I=1,NORDER SAVE=ARRAY(I,K) ARRAY(I,K)=ARRAY(I,J) 60 ARRAY(I,J)=-SAVE C ACCUMULATE ELEMENTS OF INVERSE MATRIX 61 DO 70 I=1,NORDER IF(I-K) 63,70,63 63 ARRAY(I,K)=-ARRAY(I,K)/AMAX 70 CONTINUE 71 DO 80 I=1,NORDER DO 80 J=1,NORDER IF(I-K) 74,80,74 74 IF(J-K) 75,80,75 75 ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J) 80 CONTINUE 81 DO 90 J=1,NORDER IF(J-K) 83,90,83 83 ARRAY(K,J)=ARRAY(K,J)/AMAX 90 CONTINUE ARRAY(K,K)=1./AMAX 100 DET=DET*AMAX C RESTORE ORDERING OF MATRIX 101 DO 130 L=1,NORDER K=NORDER-L+1 J=IK(K) IF(J-K) 111,111,105 105 DO 110 I=1,NORDER SAVE=ARRAY(I,K) ARRAY(I,K)=-ARRAY(I,J) 110 ARRAY(I,J)=SAVE 111 I=JK(K) IF(I-K) 130,130,113 113 DO 120 J=1,NORDER SAVE=ARRAY(K,J) ARRAY(K,J)=-ARRAY(I,J) 120 ARRAY(I,J)=SAVE 130 CONTINUE 140 RETURN END C C+ SUBROUTINE REFRACTION(ARA,ADEC,RA,DEC,EL,H) C C Calculate the true declination and right ascension from C apparent right ascension and declination by applying a C correction for atmospheric refraction. C C This routine is a modification of the Norris subroutine C rad_refract.f, which calculates apparent positions from true ones. C C The values for Mauna Kea are hard coded in. C C RA and DEC are in radians. C C Parameters - (">" input, "<" output) C C (>) ARA (Real*8) Apparent Right Ascention in radians C C (>) ADEC (Real*8) Apparent Declination in radians C C (<) RA (Real*8) True Right Ascention in radians C C (<) DEC (Real*8) True Declination in radians C C (>) H (Real*8) Hour angle in radians C C (>) LAT (Real*8) Observer's latitude in radians C C (>) PRES (Real*8) Local air pressure in mm Hg. C C (>) TEMP (Real*8) Local air temperature in Celsius. C C (>) WAVE (Real*8) Wavelength in micrometers C C (>) VP (Real*8) Vapor pressure in mm Hg (not currently used.) C C (>) EL (Real*8) Elevation in degrees C C C No external references. C C History: CRO/CIT. 20 June 1988. Original. C C Modified Oct. 1994, J.Cohen - Caltech C+ IMPLICIT NONE C REAL*8 ARA,ADEC,RA,DEC,LAT,PRES,TEMP,WAVE,VP C REAL*8 PI,RPD,A,Z,Q,R,DA,DD,HALFPI,N,NUM,TCORR,TANZ,EL,H, : WAVERS,DPR,AD,ZD,RPAS,COSQ,SINQ,SINA,SINZ,COSZ C PARAMETER (PI =3.14159265358979324, : RPD =PI/180., : HALFPI=PI/2.0, : RPAS = RPD/3600.0, : DPR = 1.0/RPD) C C Parameters for Mauna Kea C LAT=20.0*RPD WAVE=0.5 TEMP=2.0 PRES=456.0 C C Altitude and Zenith Distance C SINA = DSIN(LAT)*DSIN(ADEC) + DCOS(LAT)*DCOS(ADEC)*DCOS(H) COSZ = SINA SINZ = DSQRT(1-SINA*SINA) TANZ = SINZ/COSZ C C Parallactic angle, q C C Q = DASIN( DCOS(LAT)*DSIN(H)/DSIN(Z) ) C IF (SINZ .NE. 0.0) THEN SINQ = DCOS(LAT)*DSIN(H)/SINZ COSQ = (DCOS(HALFPI-LAT)-COSZ*DCOS(HALFPI-ADEC))/ : (SINZ*DSIN(HALFPI-ADEC)) ELSE C C Zenith distance is 0 so Q doesn't matter. R=0 anyway. C SINQ = 0.0 COSQ = 0.0 END IF C C check quadrant of q - note that correct formula for C cos(q) is expression below divided by C sin(z)*sin(ninety-dec), but only numerator C suffices to get sign of cos(q). C C cosq = (dcos(halfpi-lat) - dcos(z)*dcos(halfpi-dec)) C if (cosq.lt.0.0) then C q = pi-q c type *,' q in second quadrant,dec,ha ',dec,h C end if C C Refractive index at standard TP C WAVERS = 1./(WAVE*WAVE) N = 64.328 + (29498.1/(146.-WAVERS)) + (255.4/(41.-WAVERS)) N = 1.E-6*N C C Temperature and Pressure correction C TCORR = 1.+0.003661*TEMP NUM = 720.88*TCORR N = N*((PRES*(1.+(1.049-0.0157*TEMP)*1.E-6*PRES))/NUM) C C Correction for water vapor pressure C C N = N - (((0.0624-0.000680*WAVERS)/TCORR)*VP) C C AD = A*DPR ZD = DASIN(SINZ)*DPR C C R, refraction C R = N * 206265. * TANZ C C Changes in RA and DEC in radians. C DA = R*SINQ*RPAS/DCOS(ADEC) DD = R*COSQ*RPAS C C PRINT *,'COS(Z) SIN(Z) = ',COSZ,SINZ C PRINT *,'COS(Q) SIN(Q) = ',COSQ,SINQ C PRINT *,'DA,DD = ',DA,DD C C Corrected RA and DEC C RA = ARA - DA DEC = ADEC - DD C END C PARSING ROUTINES (higher level followed by primitive) C..higher level parsing aids C..in each case, "line, i1, i2, narg" are as in d_parse C..in each case, default a default value for null tokens C..GET_INT (parses next token and inteprets as an integer value) integer function get_int (line, i1, i2, narg, default) implicit none character*(*) line integer i1, i2, len, narg integer default integer d_parse character*80 vstring character*8 fmt narg = d_parse (line, i1, i2, vstring, len) if (len.le.0) then get_int = default return endif write (fmt, '(''(I'',I2.2,'') '')') len read (vstring(1:len), fmt, err=901) get_int return 901 write (*, '(''Error reading integer variable, near:'')') write (*, *) vstring(1:len) stop end C..GET_REAL (parses next token and inteprets as a real value) real function get_real (line, i1, i2, narg, default) implicit none character*(*) line integer i1, i2, len, narg real default integer d_parse character*80 vstring character*8 fmt narg = d_parse (line, i1, i2, vstring, len) if (len.le.0) then get_real = default return endif write (fmt, '(''(F'',I2.2,''.0) '')') len read (vstring(1:len), fmt, err=901) get_real return 901 write (*, '(''Error reading real variable, near:'')') write (*, *) vstring(1:len) stop end C..GET_TIME (parses next token and inteprets as a time/degree format) real*8 function get_time (line, i1, i2, narg, default) implicit none character*(*) line integer i1, i2, len, narg real*8 default integer d_parse character*80 vstring character*20 tstring integer tlen, ndx, tdx, colondx, j logical neg_sign integer hhh, mm real ss narg = d_parse (line, i1, i2, vstring, len) if (len.le.0) then get_time = default return endif C....rewrite vstring into tstring, using format ddd:mm:ss.s.... tstring(1:20) = ' ' ndx = 1 if (vstring(ndx:ndx).eq.'-') then neg_sign = .TRUE. vstring(ndx:ndx) = ' ' ndx = ndx + 1 else neg_sign = .FALSE. endif colondx = 0 do 100 j = ndx, len if (vstring(j:j).ne.':') go to 100 colondx = j go to 101 100 continue 101 continue if (colondx.eq.0) go to 902 tdx = 4 - (colondx - ndx) tstring(tdx:4) = vstring(ndx:colondx) ndx = colondx + 1 colondx = 0 do 200 j = ndx, len if (vstring(j:j).ne.':') go to 200 colondx = j go to 201 200 continue 201 continue if (colondx.eq.0) go to 902 tdx = 7 - (colondx - ndx) if (tdx.le.4) go to 902 tstring(tdx:7) = vstring(ndx:colondx) ndx = colondx + 1 tdx = 8 tlen = tdx + (len - ndx) tstring(tdx:tlen) = vstring(ndx:len) read (tstring(1:20), '(I3,X,I2,X,F13.0)', err=901) hhh, mm, ss get_time = hhh + mm/60. + ss/3600. if (neg_sign) get_time = -(get_time) return 901 write (*, '(''Error reading time variable, near:'')') write (*, *) vstring(1:len) stop 902 write (*, '(''Incorrect time format, near:'')') write (*, *) vstring(1:len) stop end C..GET_FIELD (parses next token(field) and places in a character string C......of size "nchar" subroutine get_field (line, i1, i2, narg, field, nchar) implicit none character*(*) line integer i1, i2, narg character*(*) field integer nchar integer d_parse character*80 vstring integer len vstring(1:20) = ' ' vstring(1:40) = ' ' vstring(1:60) = ' ' vstring(1:80) = ' ' narg = d_parse (line, i1, i2, vstring, len) if (len.lt.0) then return endif field(1:nchar) = vstring(1:nchar) return end C..D_PARSE: the parsing primitive integer function d_parse (line, i1, i2, vstring, len) C....This routine parses a whitespace-delimited line, returning one C.. argument of length "len" in the string "vstring". Strings C.. enclosed in matching quotes ("") or apostrophes (") are C.. also handled, with leading and trailing white-space removed. C.. The inputs are "line" and the search start/end characters, "i1" and "i2". C.. If "i1" is 1, the number of parsed arguments is reset. "i1" is C.. updated to the next appropriate character for further searches. C.. Note that the "#" character is interpreted to mean end-of line C.. BUGS: "whitespace" currently means blanks; tabs should be added character*(*) line integer i1, i2 character*(*) vstring integer len integer narg integer first, last, endstr integer j logical str_flag save if (i1.eq.1) narg = 0 if (i1.gt.i2) then len = -1 d_parse = narg return endif C....Search for first non-blank character first = 0 do 100 j = i1, i2 if (line(j:j).eq.' ') go to 100 first = j go to 101 100 continue 101 continue C....Check for no more args, including "#" indicating EOL (for comments) if (first.eq.0 .or. line(first:first).eq.'#') then len = -1 i1 = i2 + 1 d_parse = narg return endif C....Check for string delimiters, " and ' endstr = 0 str_flag = .FALSE. if (line(first:first).eq.'"') then first = first + 1 do 200 j = first, i2 if (line(j:j).ne.'"') go to 200 endstr = j go to 201 200 continue 201 if (endstr.eq.0) go to 901 if (endstr.lt.i2 .and. line(endstr+1:endstr+1).ne.' ') go to 902 str_flag = .TRUE. endif if (line(first:first).eq.'''') then first = first + 1 do 250 j = first, i2 if (line(j:j).ne.'''') go to 250 endstr = j go to 251 250 continue 251 if (endstr.eq.0) go to 901 if (endstr.lt.i2 .and. line(endstr+1:endstr+1).ne.' ') go to 902 str_flag = .TRUE. endif C....Check for null argument if (str_flag .and. endstr.eq.first) then i1 = endstr + 1 len = 0 narg = narg + 1 d_parse = narg return endif C....Define last character for search if (str_flag) then last = endstr - 1 else last = i2 endif C....Search for first non-blank character (redundant except for strings) do 400 j = first, last if (line(j:j).eq.' ') go to 400 first = j go to 401 400 continue 401 continue C....If string, search backward for end; otherwise, search forward if (str_flag) then do 420 j = last, first+1, -1 if (line(j:j).eq.' ') go to 420 last = j go to 500 420 continue else do 430 j = first+1, last if (line(j:j).ne.' ') go to 430 last = j - 1 go to 500 430 continue endif 500 continue C....We have a valid argument; augment counter, calculate length & copy string narg = narg + 1 len = last - first + 1 vstring(1:len) = line(first:last) d_parse = narg C....Increment new-search position; if string, skip over the second delimiter i1 = last + 1 if (str_flag) i1 = endstr + 1 return 901 write (*,'(''Unmatched string delimiters'')') stop 902 write (*,'(''End of string not followed by whitespace'')') stop end SUBROUTINE sla_DAFIN (STRING, IPTR, A, J) *+ * - - - - - - * D A F I N * - - - - - - * * Sexagesimal character string to angle (double precision) * * Given: * STRING c*(*) string containing deg, arcmin, arcsec fields * IPTR i pointer to start of decode (1st = 1) * * Returned: * IPTR i advanced past the decoded angle * A d angle in radians * J i status: 0 = OK * +1 = default, A unchanged * -1 = bad degrees ) * -2 = bad arcminutes ) (note 3) * -3 = bad arcseconds ) * * Example: * * argument before after * * STRING '-57 17 44.806 12 34 56.7' unchanged * IPTR 1 16 (points to 12...) * A ? -1.00000D0 * J ? 0 * * A further call to sla_DAFIN, without adjustment of IPTR, will * decode the second angle, 12deg 34min 56.7sec. * * Notes: * * 1) The first three "fields" in STRING are degrees, arcminutes, * arcseconds, separated by spaces or commas. The degrees field * may be signed, but not the others. The decoding is carried * out by the DFLTIN routine and is free-format. * * 2) Successive fields may be absent, defaulting to zero. For * zero status, the only combinations allowed are degrees alone, * degrees and arcminutes, and all three fields present. If all * three fields are omitted, a status of +1 is returned and A is * unchanged. In all other cases A is changed. * * 3) Range checking: * * The degrees field is not range checked. However, it is * expected to be integral unless the other two fields are absent. * * The arcminutes field is expected to be 0-59, and integral if * the arcseconds field is present. If the arcseconds field * is absent, the arcminutes is expected to be 0-59.9999... * * The arcseconds field is expected to be 0-59.9999... * * 4) Decoding continues even when a check has failed. Under these * circumstances the field takes the supplied value, defaulting * to zero, and the result A is computed and returned. * * 5) Further fields after the three expected ones are not treated * as an error. The pointer IPTR is left in the correct state * for further decoding with the present routine or with DFLTIN * etc. See the example, above. * * 6) If STRING contains hours, minutes, seconds instead of degrees * etc, or if the required units are turns (or days) instead of * radians, the result A should be multiplied as follows: * * for to obtain multiply * STRING A in A by * * d ' " radians 1 = 1D0 * d ' " turns 1/2pi = 0.1591549430918953358D0 * h m s radians 15 = 15D0 * h m s days 15/2pi = 2.3873241463784300365D0 * * Called: sla_DFLTIN * * P.T.Wallace Starlink 13 September 1990 *- IMPLICIT NONE CHARACTER*(*) STRING INTEGER IPTR DOUBLE PRECISION A INTEGER J DOUBLE PRECISION AS2R PARAMETER (AS2R=4.84813681109535993589914102358D-6) INTEGER JF,JD,JM,JS DOUBLE PRECISION DEG,ARCMIN,ARCSEC * Preset the status to OK JF=0 * Defaults DEG=0D0 ARCMIN=0D0 ARCSEC=0D0 * Decode degrees, arcminutes, arcseconds CALL sla_DFLTIN(STRING,IPTR,DEG,JD) IF (JD.GT.1) THEN JF=-1 ELSE CALL sla_DFLTIN(STRING,IPTR,ARCMIN,JM) IF (JM.LT.0.OR.JM.GT.1) THEN JF=-2 ELSE CALL sla_DFLTIN(STRING,IPTR,ARCSEC,JS) IF (JS.LT.0.OR.JS.GT.1) THEN JF=-3 * See if the combination of fields is credible ELSE IF (JD.GT.0) THEN * No degrees: arcmin, arcsec ought also to be absent IF (JM.EQ.0) THEN * Suspect arcmin JF=-2 ELSE IF (JS.EQ.0) THEN * Suspect arcsec JF=-3 ELSE * All three fields absent JF=1 END IF * Degrees present: if arcsec present so ought arcmin to be ELSE IF (JM.NE.0.AND.JS.EQ.0) THEN JF=-3 * Tests for range and integrality * Degrees ELSE IF (JM.EQ.0.AND.DINT(DEG).NE.DEG) THEN JF=-1 * Arcminutes ELSE IF ((JS.EQ.0.AND.DINT(ARCMIN).NE.ARCMIN).OR. : ARCMIN.GE.60D0) THEN JF=-2 * Arcseconds ELSE IF (JS.LT.0.OR.ARCSEC.GE.60D0) THEN JF=-3 END IF END IF END IF * Unless all three fields absent, compute angle value IF (JF.LE.0) THEN A=AS2R*(60D0*(60D0*ABS(DEG)+ARCMIN)+ARCSEC) IF (JD.LT.0) A=-A END IF * Return the status J=JF END SUBROUTINE sla_DD2TF (NDP, DAYS, SIGN, IHMSF) *+ * - - - - - - * D D 2 T F * - - - - - - * * Convert an interval in days into hours, minutes, seconds * (double precision) * * Given: * NDP int number of decimal places of seconds * DAYS dp interval in days * * Returned: * SIGN char '+' or '-' * IHMSF int(4) hours, minutes, seconds, fraction * * Notes: * * 1) NDP less than zero is interpreted as zero. * * 2) The largest useful value for NDP is determined by the size * of DAYS, the format of DOUBLE PRECISION floating-point numbers * on the target machine, and the risk of overflowing IHMSF(4). * For example, on the VAX, for DAYS up to 1D0, the available * floating-point precision corresponds roughly to NDP=12. * However, the practical limit is NDP=9, set by the capacity of * the 32-bit integer IHMSF(4). * * 3) The absolute value of DAYS may exceed 1D0. In cases where it * does not, it is up to the caller to test for and handle the * case where DAYS is very nearly 1D0 and rounds up to 24 hours, * by testing for IHMSF(1)=24 and setting IHMSF(1-4) to zero. * * P.T.Wallace Starlink 11 September 1990 *- IMPLICIT NONE INTEGER NDP DOUBLE PRECISION DAYS CHARACTER SIGN*(*) INTEGER IHMSF(4) * Days to seconds DOUBLE PRECISION D2S PARAMETER (D2S=86400D0) INTEGER NRS,N DOUBLE PRECISION RS,RM,RH,A,AH,AM,AS,AF * Handle sign IF (DAYS.GE.0D0) THEN SIGN='+' ELSE SIGN='-' END IF * Field units in terms of least significant figure NRS=1 DO N=1,NDP NRS=NRS*10 END DO RS=DBLE(NRS) RM=RS*60D0 RH=RM*60D0 * Round interval and express in smallest units required A=ANINT(RS*D2S*ABS(DAYS)) * Separate into fields AH=AINT(A/RH) A=A-AH*RH AM=AINT(A/RM) A=A-AM*RM AS=AINT(A/RS) AF=A-AS*RS * Return results IHMSF(1)=MAX(NINT(AH),0) IHMSF(2)=MAX(MIN(NINT(AM),59),0) IHMSF(3)=MAX(MIN(NINT(AS),59),0) IHMSF(4)=MAX(NINT(MIN(AF,RS-1D0)),0) END SUBROUTINE sla_DFLTIN (STRING, NSTRT, DRESLT, JFLAG) *+ * - - - - - - - * D F L T I N * - - - - - - - * * Convert free-format input into double precision floating point * * Given: * STRING c string containing number to be decoded * NSTRT i pointer to where decoding is to start * DRESLT d current value of result * * Returned: * NSTRT i advanced to next number * DRESLT d result * JFLAG i status: -1 = -OK, 0 = +OK, 1 = null, 2 = error * * Called: sla__IDCHF * * Notes: * * 1 The reason DFLTIN has separate OK status values for + * and - is to enable minus zero to be detected. This is * of crucial importance when decoding mixed-radix numbers. * For example, an angle expressed as deg, arcmin, arcsec * may have a leading minus sign but a zero degrees field. * * 2 A TAB is interpreted as a space, and lowercase characters * are interpreted as uppercase. * * 3 The basic format is the sequence of fields #^.^@#^, where * # is a sign character + or -, ^ means a string of decimal * digits, and @, which indicates an exponent, means D or E. * Various combinations of these fields can be omitted, and * embedded blanks are permissible in certain places. * * 4 Spaces: * * . Leading spaces are ignored. * * . Embedded spaces are allowed only after +, -, D or E, * and after the decomal point if the first sequence of * digits is absent. * * . Trailing spaces are ignored; the first signifies * end of decoding and subsequent ones are skipped. * * 5 Delimiters: * * . Any character other than +,-,0-9,.,D,E or space may be * used to signal the end of the number and terminate * decoding. * * . Comma is recognised by DFLTIN as a special case; it * is skipped, leaving the pointer on the next character. * See 13, below. * * 6 Both signs are optional. The default is +. * * 7 The mantissa ^.^ defaults to 1. * * 8 The exponent @#^ defaults to D0. * * 9 The strings of decimal digits may be of any length. * * 10 The decimal point is optional for whole numbers. * * 11 A "null result" occurs when the string of characters being * decoded does not begin with +,-,0-9,.,D or E, or consists * entirely of spaces. When this condition is detected, JFLAG * is set to 1 and DRESLT is left untouched. * * 12 NSTRT = 1 for the first character in the string. * * 13 On return from DFLTIN, NSTRT is set ready for the next * decode - following trailing blanks and any comma. If a * delimiter other than comma is being used, NSTRT must be * incremented before the next call to DFLTIN, otherwise * all subsequent calls will return a null result. * * 14 Errors (JFLAG=2) occur when: * * . a +, -, D or E is left unsatisfied; or * * . the decimal point is present without at least * one decimal digit before or after it; or * * . an exponent more than 100 has been presented. * * 15 When an error has been detected, NSTRT is left * pointing to the character following the last * one used before the error came to light. This * may be after the point at which a more sophisticated * program could have detected the error. For example, * DFLTIN does not detect that '1D999' is unacceptable * (on the VAX) until the entire number has been decoded. * * 16 Certain highly unlikely combinations of mantissa & * exponent can cause arithmetic faults during the * decode, in some cases despite the fact that they * together could be construed as a valid number. * * 17 Decoding is left to right, one pass. * * 18 See also FLOTIN and INTIN * * P.T.Wallace Starlink 2 August 1993 *- IMPLICIT NONE CHARACTER*(*) STRING INTEGER NSTRT DOUBLE PRECISION DRESLT INTEGER JFLAG INTEGER NPTR,MSIGN,NEXP,NDP,NVEC,NDIGIT,ISIGNX,J DOUBLE PRECISION DMANT,DIGIT * Current character NPTR=NSTRT * Set defaults: mantissa & sign, exponent & sign, decimal place count DMANT=0D0 MSIGN=1 NEXP=0 ISIGNX=1 NDP=0 * Look for sign 100 CONTINUE CALL sla__IDCHF(STRING,NPTR,NVEC,NDIGIT,DIGIT) GO TO ( 400, 100, 800, 500, 300, 200, 9110, 9100, 9110),NVEC * 0-9 SP D/E . + - , ELSE END * Negative 200 CONTINUE MSIGN=-1 * Look for first leading decimal 300 CONTINUE CALL sla__IDCHF(STRING,NPTR,NVEC,NDIGIT,DIGIT) GO TO ( 400, 300, 800, 500, 9200, 9200, 9200, 9200, 9210),NVEC * 0-9 SP D/E . + - , ELSE END * Accept leading decimals 400 CONTINUE DMANT=DMANT*1D1+DIGIT CALL sla__IDCHF(STRING,NPTR,NVEC,NDIGIT,DIGIT) GO TO ( 400, 1310, 900, 600, 1300, 1300, 1300, 1300, 1310),NVEC * 0-9 SP D/E . + - , ELSE END * Look for decimal when none preceded the point 500 CONTINUE CALL sla__IDCHF(STRING,NPTR,NVEC,NDIGIT,DIGIT) GO TO ( 700, 500, 9200, 9200, 9200, 9200, 9200, 9200, 9210),NVEC * 0-9 SP D/E . + - , ELSE END * Look for trailing decimals 600 CONTINUE CALL sla__IDCHF(STRING,NPTR,NVEC,NDIGIT,DIGIT) GO TO ( 700, 1310, 900, 1300, 1300, 1300, 1300, 1300, 1310),NVEC * 0-9 SP D/E . + - , ELSE END * Accept trailing decimals 700 CONTINUE NDP=NDP+1 DMANT=DMANT*1D1+DIGIT GO TO 600 * Exponent symbol first in field: default mantissa to 1 800 CONTINUE DMANT=1D0 * Look for sign of exponent 900 CONTINUE CALL sla__IDCHF(STRING,NPTR,NVEC,NDIGIT,DIGIT) GO TO (1200, 900, 9200, 9200, 1100, 1000, 9200, 9200, 9210),NVEC * 0-9 SP D/E . + - , ELSE END * Exponent negative 1000 CONTINUE ISIGNX=-1 * Look for first digit of exponent 1100 CONTINUE CALL sla__IDCHF(STRING,NPTR,NVEC,NDIGIT,DIGIT) GO TO (1200, 1100, 9200, 9200, 9200, 9200, 9200, 9200, 9210),NVEC * 0-9 SP D/E . + - , ELSE END * Use exponent digit 1200 CONTINUE NEXP=NEXP*10+NDIGIT IF (NEXP.GT.100) GO TO 9200 * Look for subsequent digits of exponent CALL sla__IDCHF(STRING,NPTR,NVEC,NDIGIT,DIGIT) GO TO (1200, 1310, 1310, 1310, 1310, 1310, 1310, 1310, 1300),NVEC * 0-9 SP D/E . + - , ELSE END * Combine exponent and decimal place count 1300 CONTINUE NPTR=NPTR-1 1310 CONTINUE NEXP=NEXP*ISIGNX-NDP * Skip if net exponent negative IF (NEXP.LT.0) GO TO 1500 * Positive exponent: scale up 1400 CONTINUE IF (NEXP.LT.10) GO TO 1410 DMANT=DMANT*1D10 NEXP=NEXP-10 GO TO 1400 1410 CONTINUE IF (NEXP.LT.1) GO TO 1600 DMANT=DMANT*1D1 NEXP=NEXP-1 GO TO 1410 * Negative exponent: scale down 1500 CONTINUE IF (NEXP.GT.-10) GO TO 1510 DMANT=DMANT/1D10 NEXP=NEXP+10 GO TO 1500 1510 CONTINUE IF (NEXP.GT.-1) GO TO 1600 DMANT=DMANT/1D1 NEXP=NEXP+1 GO TO 1510 * Get result & status 1600 CONTINUE J=0 IF (MSIGN.EQ.1) GO TO 1610 J=-1 DMANT=-DMANT 1610 CONTINUE DRESLT=DMANT * Skip to end of field 1620 CONTINUE CALL sla__IDCHF(STRING,NPTR,NVEC,NDIGIT,DIGIT) GO TO (1720, 1620, 1720, 1720, 1720, 1720, 9900, 1720, 9900),NVEC * 0-9 SP D/E . + - , ELSE END 1720 CONTINUE NPTR=NPTR-1 GO TO 9900 * Exits * Null field 9100 CONTINUE NPTR=NPTR-1 9110 CONTINUE J=1 GO TO 9900 * Errors 9200 CONTINUE NPTR=NPTR-1 9210 CONTINUE J=2 * Return 9900 CONTINUE NSTRT=NPTR JFLAG=J END SUBROUTINE sla_DR2AF (NDP, ANGLE, SIGN, IDMSF) *+ * - - - - - - * D R 2 A F * - - - - - - * * Convert an angle in radians to degrees, arcminutes, arcseconds * (double precision) * * Given: * NDP int number of decimal places of arcseconds * ANGLE dp angle in radians * * Returned: * SIGN char '+' or '-' * IDMSF int(4) degrees, arcminutes, arcseconds, fraction * * Notes: * * 1) NDP less than zero is interpreted as zero. * * 2) The largest useful value for NDP is determined by the size * of ANGLE, the format of DOUBLE PRECISION floating-point * numbers on the target machine, and the risk of overflowing * IDMSF(4). For example, on the VAX, for ANGLE up to 2pi, the * available floating-point precision corresponds roughly to * NDP=12. However, the practical limit is NDP=9, set by the * capacity of the 32-bit integer IDMSF(4). * * 3) The absolute value of ANGLE may exceed 2pi. In cases where it * does not, it is up to the caller to test for and handle the * case where ANGLE is very nearly 2pi and rounds up to 360 deg, * by testing for IDMSF(1)=360 and setting IDMSF(1-4) to zero. * * Called: * sla_DD2TF * * P.T.Wallace Starlink 11 September 1990 *- IMPLICIT NONE INTEGER NDP DOUBLE PRECISION ANGLE CHARACTER SIGN*(*) INTEGER IDMSF(4) * Hours to degrees * radians to turns DOUBLE PRECISION F PARAMETER (F=15D0/6.283185307179586476925287D0) * Scale then use days to h,m,s routine CALL sla_DD2TF(NDP,ANGLE*F,SIGN,IDMSF) END SUBROUTINE sla_DR2TF (NDP, ANGLE, SIGN, IHMSF) *+ * - - - - - - * D R 2 T F * - - - - - - * * Convert an angle in radians to hours, minutes, seconds * (double precision) * * Given: * NDP int number of decimal places of seconds * ANGLE dp angle in radians * * Returned: * SIGN char '+' or '-' * IHMSF int(4) hours, minutes, seconds, fraction * * Notes: * * 1) NDP less than zero is interpreted as zero. * * 2) The largest useful value for NDP is determined by the size * of ANGLE, the format of DOUBLE PRECISION floating-point * numbers on the target machine, and the risk of overflowing * IHMSF(4). For example, on the VAX, for ANGLE up to 2pi, the * available floating-point precision corresponds roughly to * NDP=12. However, the practical limit is NDP=9, set by the * capacity of the 32-bit integer IHMSF(4). * * 3) The absolute value of ANGLE may exceed 2pi. In cases where it * does not, it is up to the caller to test for and handle the * case where ANGLE is very nearly 2pi and rounds up to 24 hours, * by testing for IHMSF(1)=24 and setting IHMSF(1-4) to zero. * * Called: * sla_DD2TF * * P.T.Wallace Starlink 11 September 1990 *- IMPLICIT NONE INTEGER NDP DOUBLE PRECISION ANGLE CHARACTER SIGN*(*) INTEGER IHMSF(4) * Turns to radians DOUBLE PRECISION T2R PARAMETER (T2R=6.283185307179586476925287D0) * Scale then use days to h,m,s routine CALL sla_DD2TF(NDP,ANGLE/T2R,SIGN,IHMSF) END SUBROUTINE sla__IDCHF (STRING, NPTR, NVEC, NDIGIT, DIGIT) *+ * - - - - - - * I D C H F * - - - - - - * * Internal routine used by DFLTIN * * Identify next character in string * * Given: * STRING char string * NPTR int pointer to character to be identified * * Returned: * NPTR int incremented unless end of field * NVEC int vector for identified character * NDIGIT int 0-9 if character was a numeral * DIGIT double equivalent of NDIGIT * * NVEC takes the following values: * * 1 0-9 * 2 space or TAB !!! n.b. ASCII TAB assumed !!! * 3 D,d,E or e * 4 . * 5 + * 6 - * 7 , * 8 else * 9 outside field * * If the character is not 0-9, NDIGIT and DIGIT are either not * altered or are set to arbitrary values. * * P.T.Wallace Starlink 22 December 1992 *- IMPLICIT NONE CHARACTER*(*) STRING INTEGER NPTR,NVEC,NDIGIT DOUBLE PRECISION DIGIT CHARACTER K INTEGER NCHAR * Character/vector tables INTEGER NCREC PARAMETER (NCREC=19) CHARACTER KCTAB(NCREC) INTEGER KVTAB(NCREC) DATA KCTAB/'0','1','2','3','4','5','6','7','8','9', : ' ','D','d','E','e','.','+','-',','/ DATA KVTAB/10*1,2,4*3,4,5,6,7/ * Handle pointer outside field IF (NPTR.LT.1.OR.NPTR.GT.LEN(STRING)) THEN NVEC=9 ELSE * Not end of field: identify the character K=STRING(NPTR:NPTR) DO NCHAR=1,NCREC IF (K.EQ.KCTAB(NCHAR)) THEN * Recognized NVEC=KVTAB(NCHAR) NDIGIT=NCHAR-1 DIGIT=DBLE(NDIGIT) GO TO 2300 END IF END DO * Not recognized: check for TAB !!! n.b. ASCII assumed !!! IF (K.EQ.CHAR(9)) THEN * TAB: treat as space NVEC=2 ELSE * Unrecognized NVEC=8 END IF * Increment pointer 2300 CONTINUE NPTR=NPTR+1 END IF END SUBROUTINE sla_FLOTIN (STRING, NSTRT, RESLT, JFLAG) *+ * - - - - - - - * F L O T I N * - - - - - - - * * Convert free-format input into single precision floating point * * Given: * STRING c string containing number to be decoded * NSTRT i pointer to where decoding is to start * RESLT r current value of result * * Returned: * NSTRT i advanced to next number * RESLT r result * JFLAG i status: -1 = -OK, 0 = +OK, 1 = null, 2 = error * * Called: sla_DFLTIN * * Notes: * * 1 The reason FLOTIN has separate OK status values for + * and - is to enable minus zero to be detected. This is * of crucial importance when decoding mixed-radix numbers. * For example, an angle expressed as deg, arcmin, arcsec * may have a leading minus sign but a zero degrees field. * * 2 A TAB is interpreted as a space, and lowercase characters * are interpreted as uppercase. * * 3 The basic format is the sequence of fields #^.^@#^, where * # is a sign character + or -, ^ means a string of decimal * digits, and @, which indicates an exponent, means D or E. * Various combinations of these fields can be omitted, and * embedded blanks are permissible in certain places. * * 4 Spaces: * * . Leading spaces are ignored. * * . Embedded spaces are allowed only after +, -, D or E, * and after the decomal point if the first sequence of * digits is absent. * * . Trailing spaces are ignored; the first signifies * end of decoding and subsequent ones are skipped. * * 5 Delimiters: * * . Any character other than +,-,0-9,.,D,E or space may be * used to signal the end of the number and terminate * decoding. * * . Comma is recognised by FLOTIN as a special case; it * is skipped, leaving the pointer on the next character. * See 13, below. * * 6 Both signs are optional. The default is +. * * 7 The mantissa ^.^ defaults to 1. * * 8 The exponent @#^ defaults to E0. * * 9 The strings of decimal digits may be of any length. * * 10 The decimal point is optional for whole numbers. * * 11 A "null result" occurs when the string of characters being * decoded does not begin with +,-,0-9,.,D or E, or consists * entirely of spaces. When this condition is detected, JFLAG * is set to 1 and RESLT is left untouched. * * 12 NSTRT = 1 for the first character in the string. * * 13 On return from FLOTIN, NSTRT is set ready for the next * decode - following trailing blanks and any comma. If a * delimiter other than comma is being used, NSTRT must be * incremented before the next call to FLOTIN, otherwise * all subsequent calls will return a null result. * * 14 Errors (JFLAG=2) occur when: * * . a +, -, D or E is left unsatisfied; or * * . the decimal point is present without at least * one decimal digit before or after it; or * * . an exponent more than 100 has been presented. * * 15 When an error has been detected, NSTRT is left * pointing to the character following the last * one used before the error came to light. This * may be after the point at which a more sophisticated * program could have detected the error. For example, * FLOTIN does not detect that '1E999' is unacceptable * (on the VAX) until the entire number has been decoded. * * 16 Certain highly unlikely combinations of mantissa & * exponent can cause arithmetic faults during the * decode, in some cases despite the fact that they * together could be construed as a valid number. * * 17 Decoding is left to right, one pass. * * 18 See also DFLTIN and INTIN * * P.T.Wallace Starlink 30 March 1990 *- IMPLICIT NONE CHARACTER*(*) STRING INTEGER NSTRT REAL RESLT INTEGER JFLAG DOUBLE PRECISION DRESLT * Call the double precision version CALL sla_DFLTIN(STRING,NSTRT,DRESLT,JFLAG) IF (JFLAG.LE.0) RESLT=REAL(DRESLT) END C+ SUBROUTINE sla_INTIN (STRING, NSTRT, IRESLT, JFLAG) *+ * - - - - - - * I N T I N * - - - - - - * * Convert free-format input into an integer * * Given: * STRING c string containing number to be decoded * NSTRT i pointer to where decoding is to start * IRESLT i current value of result * * Returned: * NSTRT i advanced to next number * IRESLT i result * JFLAG i status: -1 = -OK, 0 = +OK, 1 = null, 2 = error * * Called: sla__IDCHI * * Notes: * * 1 The reason INTIN has separate OK status values for + * and - is to enable minus zero to be detected. This is * of crucial importance when decoding mixed-radix numbers. * For example, an angle expressed as deg, arcmin, arcsec * may have a leading minus sign but a zero degrees field. * * 2 A TAB is interpreted as a space. * * 3 The basic format is the sequence of fields #^, where * # is a sign character + or -, and ^ means a string of * decimal digits. * * 4 Spaces: * * . Leading spaces are ignored. * * . Spaces between the sign and the number are allowed. * * . Trailing spaces are ignored; the first signifies * end of decoding and subsequent ones are skipped. * * 5 Delimiters: * * . Any character other than +,-,0-9 or space may be * used to signal the end of the number and terminate * decoding. * * . Comma is recognised by INTIN as a special case; it * is skipped, leaving the pointer on the next character. * See 9, below. * * 6 The sign is optional. The default is +. * * 7 A "null result" occurs when the string of characters being * decoded does not begin with +,- or 0-9, or consists * entirely of spaces. When this condition is detected, JFLAG * is set to 1 and IRESLT is left untouched. * * 8 NSTRT = 1 for the first character in the string. * * 9 On return from INTIN, NSTRT is set ready for the next * decode - following trailing blanks and any comma. If a * delimiter other than comma is being used, NSTRT must be * incremented before the next call to INTIN, otherwise * all subsequent calls will return a null result. * * 10 Errors (JFLAG=2) occur when: * * . there is a + or - but no number; or * * . the number is greater than BIG (defined below). * * 11 When an error has been detected, NSTRT is left * pointing to the character following the last * one used before the error came to light. * * 12 See also FLOTIN and DFLTIN. * * P.T.Wallace Starlink 2 August 1993 *- IMPLICIT NONE CHARACTER*(*) STRING INTEGER NSTRT,IRESLT,JFLAG * Maximum allowed value DOUBLE PRECISION BIG PARAMETER (BIG=2147483647D0) INTEGER JPTR,MSIGN,NVEC,J DOUBLE PRECISION DRES,DIGIT * Current character JPTR=NSTRT * Set defaults DRES=0D0 MSIGN=1 * Look for sign 100 CONTINUE CALL sla__IDCHI(STRING,JPTR,NVEC,DIGIT) GO TO ( 400, 100, 300, 200, 9110, 9100, 9110),NVEC * 0-9 SP + - , ELSE END * Negative 200 CONTINUE MSIGN=-1 * Look for first decimal 300 CONTINUE CALL sla__IDCHI(STRING,JPTR,NVEC,DIGIT) GO TO ( 400, 300, 9200, 9200, 9200, 9200, 9210),NVEC * 0-9 SP + - , ELSE END * Accept decimals 400 CONTINUE DRES=DRES*1D1+DIGIT * Test for overflow IF (DRES.GT.BIG) GO TO 9200 * Look for subsequent decimals CALL sla__IDCHI(STRING,JPTR,NVEC,DIGIT) GO TO ( 400, 1610, 1600, 1600, 1600, 1600, 1610),NVEC * 0-9 SP + - , ELSE END * Get result & status 1600 CONTINUE JPTR=JPTR-1 1610 CONTINUE J=0 IF (MSIGN.EQ.1) GO TO 1620 J=-1 DRES=-DRES 1620 CONTINUE IRESLT=INT(DNINT(DRES)) * Skip to end of field 1630 CONTINUE CALL sla__IDCHI(STRING,JPTR,NVEC,DIGIT) GO TO (1720, 1630, 1720, 1720, 9900, 1720, 9900),NVEC * 0-9 SP + - , ELSE END 1720 CONTINUE JPTR=JPTR-1 GO TO 9900 * Exits * Null field 9100 CONTINUE JPTR=JPTR-1 9110 CONTINUE J=1 GO TO 9900 * Errors 9200 CONTINUE JPTR=JPTR-1 9210 CONTINUE J=2 * Return 9900 CONTINUE NSTRT=JPTR JFLAG=J END SUBROUTINE sla__IDCHI (STRING, NPTR, NVEC, DIGIT) *+ * - - - - - - * I D C H I * - - - - - - * * Internal routine used by INTIN * * Identify next character in string * * Given: * STRING char string * NPTR int pointer to character to be identified * * Returned: * NPTR int incremented unless end of field * NVEC int vector for identified character * DIGIT double double precision digit if 0-9 * * NVEC takes the following values: * * 1 0-9 * 2 space or TAB !!! n.b. ASCII TAB assumed !!! * 3 + * 4 - * 5 , * 6 else * 7 outside string * * If the character is not 0-9, DIGIT is either unaltered or * is set to an arbitrary value. * * P.T.Wallace Starlink 22 December 1992 *- IMPLICIT NONE CHARACTER*(*) STRING INTEGER NPTR,NVEC DOUBLE PRECISION DIGIT CHARACTER K INTEGER NCHAR * Character/vector tables INTEGER NCREC PARAMETER (NCREC=14) CHARACTER KCTAB(NCREC) INTEGER KVTAB(NCREC) DATA KCTAB/'0','1','2','3','4','5','6','7','8','9', : ' ', '+','-',','/ DATA KVTAB/10*1,2,3,4,5/ * Handle pointer outside field IF (NPTR.LT.1.OR.NPTR.GT.LEN(STRING)) THEN NVEC=7 ELSE * Not end of field: identify character K=STRING(NPTR:NPTR) DO NCHAR=1,NCREC IF (K.EQ.KCTAB(NCHAR)) THEN * Recognized NVEC=KVTAB(NCHAR) DIGIT=DBLE(NCHAR-1) GO TO 2300 END IF END DO * Not recognized: check for TAB !!! n.b. ASCII assumed !!! IF (K.EQ.CHAR(9)) THEN * TAB: treat as space NVEC=2 ELSE * Unrecognized NVEC=6 END IF * Increment pointer 2300 CONTINUE NPTR=NPTR+1 END IF END C+ subroutine dta_asfnam(string1, string2, string3, index, : string4,istatus) C C Dummy routine for LRIS coordinates program. C character*80 string1, string2, string3, string4 integer index,istatus C write(*,*) : ' You have called the dummy Figaro routine DTA_ASFNAM.' return end C+ subroutine fig_dtaerr(istatus, string1) C C Dummy routine for LRIS coordinates program. C character*80 string1 integer istatus C write(*,*)' You have called the dummy Figaro routine DTA_ERR.' return end C+ subroutine dta_rdvarf(string1,index,value,istatus) C C Dummy routine for LRIS coordinates program. C character*80 string1 integer index, istatus real*4 value C write(*,*) : ' You have called the dummy Figaro routine DTA_RDVARF.' return end C+ subroutine dta_rdvard(string1,index,dvalue,istatus) C C Dummy routine for LRIS coordinates program. C character*80 string1 integer index, istatus real*8 dvalue C write(*,*) : ' You have called the dummy Figaro routine DTA_RDVARD.' return end C+ subroutine dta_rdvarc(string1,nchar,cvalue,istatus) C C Dummy routine for LRIS coordinates program. C character*80 string1, cvalue integer nchar, istatus C write(*,*) : ' You have called the dummy Figaro routine DTA_RDVARC.' return end C+ subroutine dta_fclose(string1,istatus) C C Dummy routine for LRIS coordinates program. C character*80 string1 integer istatus C write(*,*) : ' You have called the dummy Figaro routine DTA_FCLOSE.' return end C+