C routine creates a Hydra input file using a list of x and y's plus C the header field... implicit real*8 (a-h,o-z) real*4 aweight C data files character*80 filein,xyout,radecout,radectemp character*70 out parameter(maxstars=2000) C maxstars is the maximum number of program objects that can be in the file C it is also the maximum number of possible sky positions. parameter(maxguides=2000) C "maxguides" is the maximum number of guide stars that can be in a file parameter(maxfibers=100) C "maxfibers" is the maximum number of fibers in either fiber bundle parameter(nolimit=10000) C "nolimit" is just a very big number parameter(leeslimit=4000) C leeslimit is the maximum number of objects that can be in the list C without Lee's software hic-uping. C C C C The different catagories are: objects, skys, and fops. C Each has an ra, dec, and epoch C Each has an x and y C each has an integer id number and name C in addition, the field center has an ra, dec, id, and name. C (field center is at 0,0, though) real*8 raobjs(maxstars),decobjs(maxstars),epochobjs(maxstars), . raskys(maxstars),decskys(maxstars),epochskys(maxstars), . rafops(maxguides),decfops(maxguides),epochfops(maxguides) dimension idobjs(maxstars),idskys(maxstars),idfops(maxguides) character*20 nameobjs(maxstars),nameskys(maxstars), . namefops(maxguides),namecent character*64 fieldname integer IDMSt(4),IDMSa(4) integer*4 objstatus(maxstars),skystatus(maxstars), . fopstatus(maxguides), . ixobjs(maxstars),iyobjs(maxstars), . ixfops(maxguides),iyfops(maxguides), . ixskys(maxstars),iyskys(maxstars), . itrackfops(maxguides), . itrackskys(maxstars),itrackobjs(maxstars) 1 write(*,2) 2 format(1x,'Input file name: ',$) read(*,*) filein open(unit=1,file=filein,status='old',err=111) go to 3 111 write(*,6) go to 1 3 write(*,4) 4 format(1x,'Output root name: ',$) read(*,*) out do 8 i=1,70 if(out(i:i).eq.' ') then radectemp=out(1:i-1)//'.RDTEMP' radecout=out(1:i-1)//'.RD' xyout=out(1:i-1)//'.XY' open(unit=8,file=radectemp,status='new',err=5) open(unit=9,file=radecout,status='new',err=5) open(unit=10,file=xyout,status='new',err=5) write(*,*) xyout write(8,*) xyout go to 15 5 write(*,6) write(8,6) 6 format('Im sorry, but I cant open the file. Try again') go to 3 end if 8 continue stop 'I cant handle that' C C C C Determine if this is an Old Format Nessie File C or a New Fangled Hydra File 15 continue call grabvitals(fieldname,deepoch,cuepoch,aLST,exptime, .wave,camwave,icable,aweight,1) C C C C C C C Grab the RAs and DECs of all objects, classifying them at the C same time 20 continue write(6,199) 199 format('Calling grabxy') call grabxy(racent,deccent,epochcent,idcent,namecent, . ixfops,iyfops,epochfops,idfops,namefops,Nfops,itrackfops, . maxguides, . ixskys,iyskys,epochskys,idskys,nameskys,Nskys,itrackskys, . maxstars, . ixobjs,iyobjs,epochobjs,idobjs,nameobjs,Nobjs,itrackobjs, . maxstars, . deepoch,ifile) C bug here corrected! itotal=1+Nfops+Nskys+Nobjs if(itotal.gt.leeslimit) then write(6,22) itotal,leeslimit write(8,22) itotal,leeslimit 22 format(1x,'Total number in list=',i5, 'exceeds Hydra limit',i5) stop end if C write(*,23) Nfops write(8,23) Nfops 23 format(1x,'There are',i5,' fops stars found') write(*,25) Nobjs write(8,25) Nobjs 25 format(1x,'There are',i5,' objects found') C write(*,26) Nskys write(8,26) Nskys 26 format(1x,'There are',i5,' skies found') C Now try writing them out just to see if we've got things OK: do 30 i=1, nobjs write(*,*) idobjs(i),nameobjs(i),ixobjs(i),iyobjs(i),epochobjs(i) 30 continue do 31 i=1, nskys write(*,*) idskys(i),nameskys(i),ixskys(i),iyskys(i),epochskys(i) 31 continue do 32 i=1, nfobs write(*,*) idfops(i),namefops(i),ixfops(i),iyfops(i),epochfops(i) 32 continue call dr2tf(3,racent,sign,IDMSt) call dr2af(2,deccent,sign,IDMSa) write(6,33) idcent,namecent, idmst,sign,idmsa,epochcent 33 format(i5,a20,1x,i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1, & i2.2,1x,i2.2,1x,i2.2,'.',i2.2,f7.2) C OK now we are going to do the conversions backwards... angle=0.0 write(6,123) cuepoch,aLST,exptime,camwave,wave,angle 123 format('Calling ANTI-wiynness with',6f10.2) call antiwiynness(racent,deccent,epochcent, . ixobjs,iyobjs,epochobjs,nobjs,raobjs,decobjs, . cuepoch,aLST,exptime,camwave,wave, . angle) write(*,123) cuepoch,aLST,exptime,camwave,wave,angle call antiwiynness(racent,deccent,epochcent, . ixskys,iyskys,epochskys,nskys,raskys,decskys, . cuepoch,aLST,exptime,camwave,wave, . angle) write(*,123) cuepoch,aLST,exptime,camwave,camwave,angle call antiwiynness(racent,deccent,epochcent, . ixfops,iyfops,epochfops,nfops,rafops,decfops, . cuepoch,aLST,exptime,camwave,camwave, . angle) call dr2tf(3,racent,sign,IDMSt) call dr2af(2,deccent,sign,IDMSa) write(*,101) namecent,IDMSt,sign,IDMSa,epochcent,dangle write(8,101) namecent,IDMSt,sign,IDMSa,epochcent,dangle ifred=0 write(9,1021) namecent,IDMSt,sign,IDMSa 1021 format(' 0',1x,a20,1x, .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2, .' C') Cwrite(9,101) namecent,IDMSt,sign,IDMSa,epochcent,dangle 101 format(/1x,'Field center=',a20, .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x, .i2.2,'.',i2.2, .2x,f7.2,' Ag=',f8.3) if (nobjs.lt.1) go to 110 Cwrite(9,1112) 1112 format(1x,'objects:') do 100 i=1,nobjs call dr2tf(3,raobjs(i),sign,IDMSt) call dr2af(2,decobjs(i),sign,IDMSa) write(9,1022) idobjs(i),nameobjs(i),IDMSt,sign,IDMSa write(8,102) idobjs(i),nameobjs(i),IDMSt,sign,IDMSa,epochobjs(i), . ixobjs(i),iyobjs(i),objstatus(i) 1022 format(i4,1x,a20,1x, .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2, ,' O') 102 format(i4,1x,a20, .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2, ,2x,f7.2,2i7,i3) write(10,5000) ixobjs(i),iyobjs(i) 5000 format(2i10,' O') 100 continue 110 if(nskys.lt.1) go to 1111 Cwrite(9,103) 103 format(1x,'skys:') do 1101 i=1,nskys call dr2tf(3,raskys(i),sign,IDMSt) call dr2af(2,decskys(i),sign,IDMSa) write(8,102) idskys(i),nameskys(i),IDMSt,sign,IDMSa,epochskys(i), . ixskys(i),iyskys(i),skystatus(i) write(9,1023) idskys(i),nameskys(i),IDMSt,sign,IDMSa 1023 format(i4,1x,a20,1x, .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2, ,' S') write(10,5001) ixskys(i),iyskys(i) 5001 format(2i10,' S') 1101 continue 1111 if(nfops.lt.1) go to 112 write(9,104) 104 format(1x,'fops:') do 105 i=1,nfops call dr2tf(3,rafops(i),sign,IDMSt) call dr2af(2,decfops(i),sign,IDMSa) write(8,102) idfops(i),namefops(i),IDMSt,sign,IDMSa,epochfops(i), . ixfops(i),iyfops(i),fopstatus(i) write(10,5002) ixfops(i),iyfops(i) write(8,1024) idfops(i),namefops(i),IDMSt,sign,IDMSa 1024 format(i4,1x,a20,1x, .i2.2,1x,i2.2,1x,i2.2,'.',i3.3,1x,a1,i2.2,1x,i2.2,1x,i2.2,'.',i2.2, ,' F') 5002 format(2i10,' F') 105 continue 112 continue stop 'ALL DONE!' end C C C C C This is a rewrite of "babyness" for WIYN subroutine wiynness(racent,deccent,epochcent, . ra,dec,epoch,n,ix,iy, . status, . cuepoch,aLST,exptime,camwave,wave, . dangle) implicit real*8 (a-h,o-z) real*8 ra(*),dec(*),epoch(*) integer*4 ix(*),iy(*) integer*4 status(*) dimension tempmonths(12) PARAMETER (AS2R=0.48481368110953599D-5) data encoderscale/400./ C There are 400 encoder steps per mm of travel data tempmonths/33.5,33.4,35.6,41.4,49.7,57.3,61.6,60.0,57.1, . 49.9,41.0,28.3/ month=(cuepoch-int(cuepoch))*12.+0.5 if (month.gt.12) month=month-12 if (month.lt.1) month=1 if (month.gt.12) month=12 C nighttime temperature for refraction gotten from fractional part of C current epoch tempadopt=tempmonths(month) C C C C C determine scale (encoder step per radian) C C C C call mmp(racent,deccent,0D0,0D0,epochcent,cuepoch,raz,dcz) Cwrite(9,1234) racent,deccent,epochcent 1234 format('Original plate centers are:',2g12.6,' Epoch=',f10.4) Cwrite(9,1233) raz,dcz,cuepoch 1233 format('Precessed plate centers are:',2g12.6,' Epoch=',f10.4) C get current epoch ra and dec of plate center C***************************************** C Refraction change: C iold=0 C iold=0 implies use the NEW STARLINK routines C Get refraction coefficents A and B for current wavelen call getab(tempadopt,camwave,acamwave,bcamwave,iold) call getab(tempadopt,wave,awave,bwave,iold) Cwrite(9,12322) tempadopt,camwave,acamwave,bcamwave Cwrite(9,12323) tempadopt,wave,awave,bwave 12323 format('Temp=',f5.1,' At Obj wave',f7.1,' yields refract A,B coef=', & 2g12.6) 12322 format('Temp=',f5.1,' At Cam wave',f7.1,' yields refract A,B coef=', & 2g12.6) C Now get refracted ra and dec of plate center call refractcor(raz,dcz,alst,exptime,acamwave,bcamwave,razr, .dczr,iold) C razr and dczr are now the refracted, precessed centers. C We will always use the tv wavelength for this call. C write(9,1235) razr,dczr 1235 format('Refracted plate centers are:',2g12.6) C***************************************** C C C then for each object: C C C C do 100 i=1,n C C C C C Do the precession thing Cwrite(9,9100) i, ra(i),dec(i),epoch(i) 9100 format(///'Obj',i3,' Original ra and dec=',2g12.6, & ' Epoch=',f10.4) call mmp (ra(i),dec(i),0D0,0D0,epoch(i),cuepoch,ran,dcn) C C Cwrite(9,9101) ran,dcn,cuepoch 9101 format(7x,'Precessed ra and dec=',2g12.6,' Epoch=',f10.4) C C Do the refraction thing C*************************** C Refraction change: C call refractcor(ran,dcn,alst,exptime,awave,bwave,ranr,dcnr, .iold) C Note that we are using the "wave" here; may be spectrograph or C TV depending upon subroutine call. C C C Do the geometry thing on the refracted coordinates Cwrite(9,9123) ranr,dcnr 9123 format(7x,'Refracted ra and dec=',2g12.6) call DS2TP(ranr,dcnr,razr,dczr,xi,eta) Cwrite(9,123) xi,eta 123 format('Resulting in xi, eta=',2g14.6,' (radians)') C C***************************** C C C C Correct for the corrector's pincushion and convert to mm. call correct(xi,eta,ximc,etamc) C ximc and etamc are now in mm along the curved focal plane Cwrite(9,9202) ximc,etamc 9202 format('Through the rotator (curved surface): ',2f10.3,' (mm)') C C C C Now determine where the projection of that would be on a plate plane. C These are then the x and y positions before rotation and dewarping. call flatten(ximc,etamc,ximm,etamm) Cwrite(9,9299) ximm,etamm 9299 format(5x,'Which projects to (flat surface):',2f10.3,' (mm)') C C C Correct for rotation call rotate(ximm,etamm,x,y,dangle) Cwrite(9,9203) x,y,dangle 9203 format(27x,'Rotated to:',2f10.3,' for angle=',f8.3) C Figure out where we should plunk it down on the UNWARPED board! call dewarp(x,y,xp,yp) Cwrite(9,9204) xp,yp 9204 format(22x,'Unwarped plate =',2f10.3) C convert to encoder steps xpp=encoderscale*xp ypp=encoderscale*yp C wow! Now that wasn't so bad, was it! C C C write(9,9205) xpp,ypp C9205 format(10x,'Encoder steps=',2f10.0/) C C ix(i)=xpp+0.5 iy(i)=ypp+0.5 100 continue return end C C subroutine antiwiynness(racent,deccent,epochcent, . ix,iy,epoch,n,ra,dec, . cuepoch,aLST,exptime,camwave,wave, . dangle) implicit real*8 (a-h,o-z) real*8 ra(*),dec(*),epoch(*) integer*4 ix(*),iy(*) dimension tempmonths(12) PARAMETER (AS2R=0.48481368110953599D-5) data encoderscale/400./ C There are 400 encoder steps per mm of travel data tempmonths/33.5,33.4,35.6,41.4,49.7,57.3,61.6,60.0,57.1, . 49.9,41.0,28.3/ C First, find the temperature for refraction month=(cuepoch-int(cuepoch))*12.+0.5 if (month.gt.12) month=month-12 if (month.lt.1) month=1 if (month.gt.12) month=12 C nighttime temperature for refraction gotten from fractional part of C current epoch tempadopt=tempmonths(month) C Precess plate center to current epoch: C C call mmp(racent,deccent,0D0,0D0,epochcent,cuepoch,raz,dcz) C C C write(8,1234) racent,deccent,epochcent 1234 format('Original plate centers are:',2g12.6,' Epoch=',f10.4) write(8,1233) raz,dcz,cuepoch 1233 format('Precessed plate centers are:',2g12.6,' Epoch=',f10.4) C get current epoch ra and dec of plate center C Refraction change: C iold=0 C iold=0 implies use the NEW STARLINK routines C Get refraction coefficents A and B for current wavelen call getab(tempadopt,camwave,acamwave,bcamwave,iold) call getab(tempadopt,wave,awave,bwave,iold) write(8,12322) tempadopt,camwave,acamwave,bcamwave write(8,12323) tempadopt,wave,awave,bwave 12323 format('Temp=',f5.1,' At Obj wave',f7.1,' yields refract A,B coef=', & 2g12.6) 12322 format('Temp=',f5.1,' At Cam wave',f7.1,' yields refract A,B coef=', & 2g12.6) C Now get refracted ra and dec of plate center call refractcor(raz,dcz,alst,exptime,acamwave,bcamwave,razr, .dczr,iold) C razr and dczr are now the refracted, precessed centers. C We will always use the tv wavelength for this call. write(8,1235) razr,dczr 1235 format('Refracted plate centers are:',2g12.6) C***************************************** C C C then for each object: C C C C do 100 i=1,n C C C Convert to mm in floating point: xp=ix(i)/encoderscale yp=iy(i)/encoderscale write(8,123) ix(i),iy(i),xp,yp 123 format('ANTI-wiynness:'/' Convesion to mm:',2i10,' --->', & 2f12.4) C C Next, go through Sam's dewarp routine...backwards call antidewarp(xp,yp,x,y) write(8,124) xp,yp,x,y 124 format(' Conversion through dewarp:',2f12.4,'-->',2f12.4) C Now derotate these x and y's amd give me ximm and etamm call antirotate(x,y,ximm,etamm,dangle) write(8,125) x,y,ximm,etamm,dangle 125 format('Conversion through angle:',2f12.4,'--->',2f12.4, & ' ang=',f10.2) call antiflatten(ximm,etamm,ximc,yimc) write(8,126) ximm, etamm,ximc,yimc 126 format('Conversion through flatten:',2f12.4,'--->',2f12.4) C Correct for pincushion and plate scale: call anticorrect(ximc,yimc,xi,eta) write(8,127) ximc,yimc,xi,eta 127 format('Conversion through correct:',2f16.8,2g16.8) C Go from Standard Coordinates to RA, DEC: call DTP2S(xi,eta,razr,dczr,ranr,dcnr) write(8,128) xi,eta,razr,dczr,ranr,dcnr 128 format('Converting from STD COORDS',2g16.8,' using CNT',2g16.8, & '--->',2g16.8) call antirefract(ranr,dcnr,alst,exptime,awave,bwave,ran,dcn) write(8,129) ranr,dcnr,alst,exptime,awave,bwave,ran,dcn 129 format('Correcting for refraction: coords:',2g16.8/ & 'Using LST, exptime,A, and B:',2f10.2,2g12.6/ & 'Results in ra and dec of',2g16.8) C Finally, let's precess; call mmp(ran,dcn,0D0,0D0,cuepoch,epoch(i),ra(i),dec(i)) write(8,130) ran,dcn,cuepoch,ra(i),dec(i),epoch(i) 130 format('Precession from',2g16.8,f10.2/ & 'Resuls in:',2g16.8) 100 continue return end C C C C subroutine antirotate(x,y,ximm,etamm,dangle) implicit real*8 (a-h,o-z) C note that this angle is in degrees. C CONVENTIONS C remember that +eta is +N C remember that xi is +E C we also know that when the angle is 0, fiber 216 (0,-big) should be C north C so, when angle is 0, +eta (N) should be at x= 0, y=-eta*scale -----> C C Using alegbra.... etamm=x*dcosd(dangle-90.)+y*dsind(dangle-90.) ximm=x*dsind(dangle-90.)-y*dcosd(dangle-90.) return end C subroutine rotate(ximm,etamm,x,y,dangle) implicit real*8 (a-h,o-z) C note that this angle is in degrees. C CONVENTIONS C remember that +eta is +N C remember that xi is +E C we also know that when the angle is 0, fiber 216 (0,-big) should be C north C so, when angle is 0, +eta (N) should be at x= 0, y=-eta*scale -----> C x=etamm*dcosd(dangle-90.)+ximm*dsind(dangle-90.) y=-ximm*dcosd(dangle-90.)+etamm*dsind(dangle-90.) return end C C subroutine for getting the ra and dec out of the file C C subroutine grabxy(racent,deccent,epochcent,idcent,namecent, . ixfops,iyfops,epochfops,idfops,namefops,Nfops,itrackfops, . maxfops, . ixskys,iyskys,epochskys,idskys,nameskys,Nskys,itrackskys, . maxskys, . ixobjs,iyobjs,epochobjs,idobjs,nameobjs,Nobjs,itrackobjs, . maxobjs, . deepoch,ifile) implicit real*8 (a-h,o-z) real*8 epochobjs(*), . epochskys(*), . epochfops(*) integer*4 idobjs(*),idskys(*),idfops(*), . itrackfops(*),itrackskys(*), . itrackobjs(*), . ixobjs(*),iyobjs(*), . ixskys(*),iyskys(*), . ixfops(*),iyfops(*) character*20 nameobjs(*),nameskys(*), . namefops(*),namecent character*80 header character*1 sign character*1 class C C C C C C initialize counters nobjs=0 nfops=0 nskys=0 icount=-1 1000 rewind(1) C major changes made here July 18, 1991 to allow for "comment" field and C for header field. The following restrictions apply: C a "header field" has a ":" somewhere in the first 17 columns C a comment field has a # in the first column. C everything else is assume to conform to the correct format 1001 read(1,'(A)',end=3003) header Cwrite(*,*) header if(header(1:1).eq.'#') then C comment! if(header(2:2).eq.'#') go to 1001 C keep the comment only if there is ONE #; ignore if there are two. write(9,5001) header write(10,5001) header 5001 format(a80) go to 1001 end if do 33 i=1,17 if(header(i:i).eq.':') go to 1001 33 continue C it's a comment if there is a ":" in the first 17 characters C C C OK, it's time to get serious. First find the center. Needs to be the first C object. if(header.eq.'') go to 1001 read(header(22:22),*,end=3003) class Cwrite(*,*) class if((class.eq.'C').or.(class.eq.'c')) then read(header,1003) ihr,imn,sec,sign,ideg,imin,arcsec,epochcent 1003 format(34x,i2,1x,i2,1x,f6.3,1x,a1,i2,1x,i2,1x,f5.2,2x,f7.2) racent=STOR(ihr,imn,sec) deccent=ASTOR(sign,ideg,imin,arcsec) idcent=0 namecent=header(15:34) ncent=1 go to 1001 end if C Now read in each object read(header,1002) ix,iy,class 1002 format(2I10,1x,a1) icount=icount+1 C count the number of objects, staring with 0 as the first one if((class.eq.' ').or.(class.eq.'').or. & (class.eq.'o').or.(class.eq.'O')) then C it's an object! nobjs=nobjs+1 if(nobjs.gt.maxobjs) stop 'too many objects' ixobjs(nobjs)=ix iyobjs(nobjs)=iy epochobjs(nobjs)=deepoch idobjs(nobjs)=icount write(nameobjs(nobjs),1004) icount 1004 format('FakeObj-',i4.4) itrackobjs(nobjs)=icount ELSE if ((class.eq.'F') .or.(class.eq.'f'))then nfops=nfops+1 If (nfops.gt.maxfops) stop 'too many fops' ixfops(nfops)=ix iyfops(nfops)=iy epochfops(nfops)=deepoch idfops(nfops)=icount write(namefops(nfops),1005) icount 1005 format('FakeFobs-',i4.4) itrackfops(nfops)=icount ELSE if((class.eq.'S').or.(class.eq.'s')) then nskys=nskys+1 If (nskys.gt.maxskys) stop 'too many skies' ixskys(nskys)=ix iyskys(nskys)=iy epochskys(nskys)=deepoch idskys(nskys)=icount write(nameskys(nskys),1006) icount 1006 format('FakeSkys-',i4.4) itrackskys(nskys)=icount ELSE END IF go to 1001 3003 Continue if(ncent.lt.1) stop 'cannot find the center' return end C C C C FUNCTION STOR(IHR,IMN,SEC) implicit real*8 (a-h,o-z) parameter(s2r=0.72722052166430398D-4) STOR=s2r*(60D0*(60D0*IHR+IMN)+SEC) return end FUNCTION ASTOR(SIGN,IDEG,IMIN,ARCSEC) implicit real*8 (a-h,o-z) character*1 SIGN parameter(as2r=0.48481368110953599D-5) if(SIGN.eq.'-') then ideg=-ideg imin=-imin arcsec=-arcsec end if ASTOR=as2r*(60d0*(60d0*ideg+imin)+arcsec) return end subroutine grabvitals(fieldname,deepoch,cuepoch,aLST,explength, .wave,camwav,icable,aweight,ifile) C get the field name, default epoch, current epoch, LST, exp length, C wavelength, and camera wavelength implicit real*8 (a-h,o-z) real*4 aweight C setup the names used in the hydra file character*16 Things(9),ident character*4 color(3) character*80 X character*64 fieldname,tempname logical iok(7) character*4 cable character*4 weight character*6 wts(3) C cheat cheat cheat cheat, to make reading in the data easy! C in the 'C' version of course I am sure that it will be far more C clever than this...... data Things /'FIELD NAME: ', . 'INPUT EPOCH: ', . 'CURRENT EPOCH: ', . 'SIDEREAL TIME: ' , . 'EXPOSURE LENGTH:' , . 'WAVELENGTH: ', . 'CABLE: ' , . 'GUIDEWAVE: ' , . 'WEIGHTING: ' / data color/'RED ','BLUE','PINK'/ data wts/'NONE ','WEAK ','STRONG'/ C global defaults: camwav=6000. icable=1 aweight=1.0 rewind(1) if(ifile.eq.0) then read(1,1,end=999,err=999) . fieldname(1:10),deepoch,cuepoch,aLST,explength,wave 1 format(a10,5f8.2) C format as per Nessie manual do 13 i=11,64 fieldname(i:i)=' ' 13 continue go to 100 end if C Read Hydra file 11 read(1, '(a)' ,end=99) X read(X, '(a16)' ) ident C If it's the old format, it's OK still. if(ident(1:14).eq.'SIDERIAL TIME:') .ident(1:14)='SIDEREAL TIME:' if(ident(1:11).eq.Things(1)(1:11)) then read(X, '(11X,a64)') fieldname if(fieldname(1:1).eq.' ') then tempname=fieldname(2:64) fieldname=tempname end if iok(1)=1 end if if(ident(1:12).eq.Things(2)(1:12)) then read(X, '(12X,f12.2)') deepoch if((deepoch.lt.1855.).or.(deepoch.gt.2000.)) stop . 'trouble with input epoch' iok(2)=1 end if if(ident(1:14).eq.Things(3)(1:14)) then C Big Trouble if colon here! read(X, '(14X,f12.2)') cuepoch if((cuepoch.lt.1995.).or.(cuepoch.gt.2010.)) stop . 'trouble with current epoch' iok(3)=1 end if if(ident(1:14).eq.Things(4)(1:14)) then read(X, '(14X,f12.2)') alst if((alst.lt.0.).or.(alst.gt.23.99999)) stop . 'trouble with LST entry' iok(4)=1 end if if(ident.eq.Things(5)) then read(X, '(16X,f12.2)') explength if((explength.lt.0.).or.(explength.gt.10.)) stop . 'trouble with exposure length' iok(5)=1 end if if(ident(1:11).eq.Things(6)(1:11)) then read(X, '(11X,f12.2)') wave if(wave.le.1.) wave=wave*10000. if((wave.gt.10000.).or.(wave.lt.3200.)) stop . 'trouble with wavelength' iok(6)=1 end if if(ident(1:10).eq.Things(8)(1:10)) then read(X,'(10X,f12.2)') camwav if(camwav.le.1.) camwav=camwav*10000. if((camwav.gt.10000.).or.(camwav.lt.3500.)) stop . 'trouble with guide wavelength' end if if(ident(1:6).eq.Things(7)(1:6)) then read(X,'(6X,a4)') cable if(cable(1:1).eq.' ') then if((cable(2:4).ne.'RED').and.(cable(2:4).ne.'BLU').and. . (cable(2:4).ne.'PIN')) stop . 'trouble with fiber bundle color' if(cable(2:4).eq.'BLU') icable=2 if(cable(2:4).eq.'PIN') icable=3 go to 333 end if if((cable(1:3).ne.'RED').and.(cable(1:4).ne.'BLUE').and. . (cable(1:4).ne.'PINK'))stop . 'trouble with fiber bundle color' if(cable.eq.'BLUE') icable=2 if(cable.eq.'PINK') icable=3 333 iok(7)=1 end if if(ident(1:10).eq.Things(9)(1:10)) then read(X,'(10X,a4)') weight if(weight(1:1).eq. ' ') then if((weight.ne.' STR').and.(weight(2:4).ne.'WEA').and. . (weight(2:4).ne.'NON')) stop 'dont understand this weight' if(weight.eq.' STR') aweight=1.0 if(weight(2:4).eq.'NON') aweight=0.0 if(weight(2:4).eq.'WEA') aweight=0.5 go to 11 end if if((weight.ne.'STRO').and.(weight.ne.'WEAK').and. . (weight.ne.'NONE')) stop 'dont understand this weight' if(weight.eq.'STRO') aweight=1.0 if(weight.eq.'NONE') aweight=0.0 if(weight.eq.'WEAK') aweight=0.5 end if go to 11 99 do 1101 i=1,7 if(iok(i).ne.1) then write(*,98) Things(i) 98 format(1x,'In new file could not find ',a16) stop end if 1101 continue 100 write(9,101) Things(1)(1:11),fieldname, . Things(2)(1:12),deepoch,Things(3)(1:14),cuepoch, . Things(4)(1:14),alst, . Things(5),explength,Things(6)(1:11),wave, . Things(7)(1:6),color(icable),Things(8)(1:10),camwav, . Things(9),wts(int(2.1*aweight+1)) write(10,101) Things(1)(1:11),fieldname, . Things(2)(1:12),deepoch,Things(3)(1:14),cuepoch, . Things(4)(1:14),alst, . Things(5),explength,Things(6)(1:11),wave, . Things(7)(1:6),color(icable),Things(8)(1:10),camwav, . Things(9),wts(int(2.1*aweight+1)) 101 format(a11,1x,a64,/a12,1x,f7.2,/a14,1x,f7.2,/a14,1x,f5.2,/a16,1x, . f5.2,/ . a11,1x,f5.0,/a6,1x,a4,/a10,1x,f5.0,/a10,1x,a6) write(6,1000) 1000 format('Returning from grabvitals') return 999 stop 'problem reading Nessie-format input file; I quit!' end C C C C C C C subroutine filetype(ifile) C ifile=0 Nessie file C ifile=1 Hydra file character*16 Thing ,ident data Thing /'FIELD NAME: '/ C Lee's file specifications claim that these two key words must be there C in a hydra file, but that spaces are ignored and capitization is ignored; C neither one of those things is true here, however! C we have also cheated, in so far as we are always writing these C keywords out as a16 and they start in column 1. ifile=0 rewind(1) C assume it is a Nessie file, unless proven otherwise! 1 read(1,2,end=99) ident 2 format(a16) if(ident(1:11).eq.Thing(1:11)) then ifile=1 return end if go to 1 99 return end C C C C********************************************** C Refraction change: refraction subroutine subroutine getab(temp,w,a,b,iold) real*8 colorfact(14) implicit real*8 (a-h,o-z) data alat/31.9633/ parameter(as2r=0.48481368110953599D-5) DATA EPS/1d-10/ C precision to terminate iteration (radian) DATA AS/0.0065D0/ C temperature lapse rate in the troposphere (degK/meter) DATA UPS/.6/ C relative humidity DATA h0/2064./ C height of Kitt Peak DATA p0/795./ C typical pressure on Kitt Peak in millibars data colorfact/1.046,1.026,1.014,1.005,1.000,0.996, .0.993,0.991,0.989,0.987,0.986,0.985,0.983,0.983/ c multiply refraction at 5000A by these amounts for refraction c at 3000,3500,4000,4500,5000,5500,6000,6500,7000,7500,8000,8500,9000 A c factor by which atmospheric pressure is reduced; i.e., c atmospheric pressure is 760 * P mm Hg C check on OK wavelength w: if(w.lt.3000) w=3000. if(w.gt.10000) w=10000. C then: if(iold.eq.1) then C determine coeficients the same way as before: Atemp=58.3 Btemp=-.067 C find wavelength index: i=(w-3000.)/500.+1.5 PTERM=exp(-h0/8000.) RTERM=colorfact(i)*PTERM/(0.962+0.0021*(TEMP-32.)) A=Atemp*Rterm*as2r B=Btemp*Rterm*as2r else C determine coeficients using the new & improved routines t0=273.+(temp-32.)*100./(212.-32.) C turn wavelength into microns wl=w/10000. C turn latitude into radians realat=alat*3600.*as2r C write(9,1235) h0,t0,p0,ups,wl,realat,as,eps 1235 format(1x,'Calling STARLINK with',8g12.6) call sla_REFCO(h0,t0,p0,ups,wl,realat,as,eps,a,b) end if C a and b are in radians C write(9,1234) w,a,b 1234 format(1x,'refraction coeficients for wavelength:',f10.2,2g12.6) return end subroutine antirefract(raout,decout,alst,exp,a,b,rain,decin) C Lets just do it the silly way... implicit real*8 (a-h,o-z) inum=0 rain=raout decin=decout write(8,*) raout,decout,alst,exp,a,b 12 call refractcor(rain,decin,alst,exp,a,b,ratest,dectest,0) inum=inum+1 difra=(ratest-raout)*206265. difdec=(dectest-decout)*206265. write(8,13) inum,difra,difdec 13 format('Antirefract:',i5,2f10.2) if((abs(difra).lt.0.01).and.(abs(difdec).lt.0.01)) return rain=rain-difra/2./206265. decin=decin-difdec/2./206265. go to 12 end C C C C C C This subroutine replaces "aver" subroutine refractcor(ra,dec,st,exptime,a,b,rar,decr,iold) implicit real*8 (a-h,o-z) real*8 r(10),d(10) Cwrite(9,999) ra,dec,st,exptime,a,b 999 format(1x,'Doing refractin correction for ra,dec,lst,exp,a,b=', ./1x,6g12.6) zo=dacosd(dsin(dec)*dsind(31.9633d0)+dcos(dec)*dcosd(31.9633d0)* . dcos(st*15.*1.745329d-2-ra)) division=exptime/10. sts=st-(4.5*division) SUMRA=0. SUMDC=0. SUMWEIGHTS=0. do 1 i=1,10 stnow=sts+(i-1)*division call refract(ra,dec,stnow,a,b,rra,rdc,AIR,iold) r(i)=rra d(i)=rdc SUMRA=SUMRA+r(i)/AIR SUMDC=SUMDC+d(i)/AIR SUMWEIGHTS=SUMWEIGHTS+1./AIR 1 CONTINUE rave=sumra/sumweights dcave=sumdc/sumweights C OK, now correct the RA's and Dec's rar=ra+rave decr=dec+dcave zp=dacosd(dsin(decr)*dsind(31.9633d0)+dcos(decr)*dcosd(31.9633d0)* . dcos(st*15.*1.745329d-2-rar)) dz=(zo-zp)*3600. Cwrite(9,123) ra,rar,dec,decr,rave,dcave,zo,zp,dz C23 format(1x,'REFRACTION ra,dec=', C .1x,4g12.6,/1x,2g12.4,2f10.4,f6.1) return end subroutine refract(ra,dec,st,a,b,rra,rdc,air,iold) c INPUT: c ra is ra in radians c dec is dec in radians c st is siderial time in decimal hrs C a and b are the refraction coefficients for r=A*tanz + Btan^3z c OUTPUT: c rra is correction for refraction in right ascention in RADIANS c corrected for declination c rdc is correction for refraction in declination in RADIANS c air is the airmass implicit real*8 (a-h,o-z) DATA degtorad/0.0174532925d0/ data pi/3.141592653589793d0/ data alat/31.9633/ realat=alat*pi/180. twopi=pi*2. RST=ST*15.*pi/180. HAold=RST-RA if(haold.ge.pi)haold=haold-twopi if(haold.lt.-pi)haold=haold+twopi C see to it that ha is between -pi and +pi cz=dsin(dec)*dsin(realat)+dcos(dec)*dcos(realat)*dcos(HAold) air=1./cz z=dacos(cz) Cwrite(9,991) ra,dec,st,a,b,z,air C91 format(1x,'unrefracted ra,dec,st,a,b,z,air=', C .6f12.4,f10.2) C zold is now the old zenith distance if(iold.eq.1) then zr=z-A*dtan(z)+B*(dtan(z))**3. C this is almost right. A better thing would be to do dz=(z-zr)*206280.d0 Cwrite(9,7191) z,zr,dz C191 format(1x,'With OLDCODE: z, zr, dz=',2g12.6,f10.2) else call sla_refz(z,a,b,zr) dz=(z-zr)*206280.d0 Cwrite(9,7192) z,zr,dz C192 format(1x,'With STARLINK: z, zr, dz=',2g12.6,f10.2) end if C znew is now the new zenith distance Cwrite(9,990) zr C90 format(1x,'Leading to new zenith distance=',f12.4) x=dsin(HAold)/(dcos(HAold)*dsin(realat)-dtan(dec)*dcos(realat)) C x is now the arctan of the azimuth: if(Haold.ge.0.) then if(x.ge.0.0d0) then Az=pi+datan(x) else Az=twopi+datan(x) end if else if(x.ge.0.d0) then Az=datan(x) else Az=pi+datan(x) end if end if 811 format(1x,'z,zr,HAold,decold,decnew,x,Az=', ./7g12.6) C Az should be between 0 and twopi at this point. C find new declination: s=dsin(realat)*dcos(zr)+dcos(realat)*dsin(zr)*dcos(Az) decnew=dasin(s) Cwrite(9,811)z,zr,Haold,dec,decnew,x,Az rdc=decnew-dec C find new hour hangle: sinha=-dsin(Az)*dsin(zr)/dcos(decnew) cosha=(dcos(zr)-dsin(decnew)*dsin(realat))/ . (dcos(decnew)*dcos(realat)) hanew=datan2(sinha,cosha) rra=haold-hanew if(rra.gt.pi) rra=rra-twopi if(rra.lt.-pi) rra=rra+twopi Cwrite(9,780) rra,rdc C80 format(1x,'Refraction=',2g12.6) C change has to be between -pi and pi (and of course is usually much smaller!) return end subroutine anticorrect(ximc,etamc,xi,eta) implicit real*8 (a-h,o-z) * ARC SECONDS TO RADIANS PARAMETER (AS2R=4.8481368110953599D-6) Cdata scale/9.4032/ data scale/9.4000/ C xout is in mm xout=ximc yout=etamc C xi and eta are in radians xi=xout*scale*as2r eta=yout*scale*as2r C Start with the assumption that the output is similar to the input C xi and eta are in radians; xtest and ytest are in mm 12 call correct(xi,eta,xtest,ytest) write(6,*) xi,eta,xtest,ytest inum=inum+1 C difx is in mm difx=xtest-xout dify=ytest-yout if ((abs(difx).lt. 0.00001).and.(abs(dify).lt. 0.00001)) return C Take half the difference C Put it back in radians: xi=xi-(difx/2.)*scale*as2r eta=eta-(dify/2.)*scale*as2r go to 12 end subroutine correct(xi,eta,ximc,etamc) C this program considerably simplied after our initial astrometry 3/95 C old code: C implicit real*8 (a-h,o-z) C data AA,BB,CC,DD/6.391006d0,1.249909d-4,2.576783d-5,4.136236d-7/ C d=dsqrt(xi**2+eta**2) C dp=AA*d+BB*d**2+CC*d**3+DD*d**4 C Great. Now, what are the new xi's and eta's? C ximc=xi*(dp/d) C etamc=eta*(dp/d) C ximc and etamc are now xi and eta in the unlikely units of "mm" C implicit real*8 (a-h,o-z) * ARC SECONDS TO RADIANS PARAMETER (AS2R=4.8481368110953599D-6) Cdata scale/9.4032/ data scale/9.4000/ C this done for Sam 8/28/95 data c/81.0/ C Correct for pincushion distortion: f=1.+c*(xi**2+eta**2) xip=xi*f etap=eta*f Cwrite(9,101) f,xip,etap 101 format('Pincushion factor=',f10.6/ & 'Resulting in xi, eta=',2g14.6) C Correct for plate scale ximc=xip/scale/as2r etamc=etap/scale/as2r C OK: ximc and etamc are now eta and xi in mm along curved surface return end C C subroutine antidewarp(xout,yout,xin,yin) implicit real*8 (a-h,o-z) C Cheat! Just keep calling dewarp until we get good answers: inum=0 xin=xout yin=yout C Start with the assumption that xin and yin are pretty close to C to the final answers: 12 call dewarp(xin,yin,xtest,ytest) inum=inum+1 difx=xtest-xout dify=ytest-yout write(8,13) inum,difx,dify 13 format('antidewarp:',i5,2f10.3) if ((abs(difx).lt.0.0001).and.(abs(dify).lt.0.0001)) return C Take half the difference xin=xin-difx/2. yin=yin-dify/2. go to 12 end C warp subroutine: subroutine dewarp(xin,yin,xout,yout) implicit real*8 (a-h,o-z) C These are Sam's estimage of the radius from the pupil to the focal C surface data rfocal/5675.2/ C C What is the actual orthogonal distance that we want from the C radius vector to the top of the fiber. din=dsqrt(xin**2+yin**2) C The following were modified as per Sam's 13 March note C re the "mid radius" where length is conserved.... if(din.lt.151.9325) then reff=5679.881 else temp=din*400. reff=5670.799 + temp*(2.988726e-4 - 2.458933e-9 * temp) C Note that Sam's equation was in terms of encoder units. end if C Great. Now, what is the distance we want along the plate prior to C being warped? dout=reff*dasin(din/rfocal) C Now make sure that the proportions work out: xout=0. yout=0. if(din.gt.0.) then yout=(dout/din)*yin xout=(dout/din)*xin end if Cwrite(9,1234) din,dout,reff,rfocal,xin,xout,yin,yout 1234 format('Din='f10.3,' Dout=',f10.3,' reff=',f10.3, 'rfoc=', & f10.3,'xin,out; yin,out=',f10.3) return end C Anti-deflatten! subroutine antiflatten(xout,yout,xin,yin) implicit real*8 (a-h,o-z) xin=xout yin=yout C Start with the assumption that the output is similar to the input 12 call flatten(xin,yin,xtest,ytest) inum=inum+1 difx=xtest-xout dify=ytest-yout write(8,13) inum,difx,dify 13 format('antiflatten:',i5,2f10.3) if ((abs(difx).lt.0.0001).and.(abs(dify).lt.0.0001)) return C Take half the difference xin=xin-difx/2. yin=yin-dify/2. go to 12 end C C This subroutine re-written and considerably simplied from old Hydra C days! subroutine flatten(x,y,xp,yp) implicit real*8 (a-h,o-z) r=5675.2 C = 223.4 inches radius of curvature c r is radius of curve (mm) of Lockheed camera plate adapted for WIYN zed=sqrt(x**2+y**2) theta=zed/r h=r*sin(theta) xp=0. yp=0. if (zed.gt.0.) then factor=h/zed xp=x*factor yp=y*factor end if return end C C subroutines used (from kaboomlib.a) C C call libinit() [no arguments]: initializes the "data base"; use before any C other of the C calls C C call rotmode (flag) turns rotation on (1) or off (0). C C call rotangle (real*8 angle in radians): sets the rotation angle C C call concenmode (flag) turns concentricity correction on (1) or off (0) C C call setfib (ibut, x,y,istat) sets button ibut to integer "sky" coordinates C x,y. Returns istat (integer) 0 = sucessful; 1 = if fiber would be over- C extended or pushed back beyond the park circle; 2= if the attempted move was C to a step C C call checkfib(ibut,x,y,istat) checks to see if button ibut can go to position C x,y; same as setfib but the database is not updated. C C call setfibphys/checkfibphys [same arguments as above]: does the same but C uses PHYSICAL rather than SKY coordinates. I use this if whereis has C returns the x and y position of a button that I may need to put back C there. C C call whereis (ibut,x1,y1,x2,y2,x3,y3) : returns the location (x1,y1) of C button ibut [ PHYSCIAL coordinates]. x2,y2 should always be the same C as (x1,y1), and are the coordinates used for collision testing C (x3,y3) are the taget location for the next move and are only sometimes C valid C C call park (ibut,flag) Parks the fiber; 0 = park; 1 = stow position C C ouch=interract (ibuta,ibutb) 0 = they aren't in each other's personal C space; 1=they are C C ouch=ibuttonhitsanother returns 0 if there is no interraction with C another button; !=0 if there is. I don't use this. C C the rest of these stuff is directly from the STARLINK library SUBROUTINE sla_REFRO (Z0,H0,T0,P0,UPS,WL,PH,AS,EPS,REF) *+ * * - - - - - - * R E F R O * - - - - - - * * * Atmospheric refraction for optical wavelengths * * * Given: * Z0 dp observed zenith distance of the source (radian) * H0 dp height of the observer above sea level (metre) * T0 dp ambient temperature at the observer (deg K) * P0 dp pressure at the observer (millibar) * UPS dp relative humidity at the observer (range 0-1) * WL dp effective wavelength of the source (micrometre) * PH dp latitude of the observer (radian, astronomical) * AS dp temperature lapse rate in the troposphere (degK/metre) * EPS dp precision required to terminate iteration (radian) * * Returned: * REF dp refraction: in vacuo ZD minus observed ZD (radian) * * Called: sla_$ATMT, sla_$ATMS * * Typical values for the AS and EPS arguments would be 0.0065D0 and * 1D-10 respectively. * * This routine computes the refraction for zenith distances up to * and beyond 90 deg using the method of Hohenkerk and Sinclair (NAO * Technical Notes 59 and 63). The code is a trivially modified form * of the AREF subroutine of C.Hohenkerk (HMNAO, September 1984). * Most of the modifications are cosmetic; in addition the angle * arguments have been changed to radians and the Z0 = 0 case is * handled. For Z0 > 90 deg it is up to the caller to supply an * H0 consistent with the depression of the horizon thus implied. * * Fixed values of the water vapour exponent, height of tropopause, and * height at which refraction is negligible are used. * * Original version by C.Hohenkerk, HMNAO, September 1986. * Starlink version by P.T.Wallace, July 1986. * *+ IMPLICIT NONE DOUBLE PRECISION Z0,H0,T0,P0,UPS,WL,PH,AS,EPS,REF * * Fixed parameters * DOUBLE PRECISION GCR,DMD,DMW,S,GAMMA,HT,HS * Universal gas constant PARAMETER (GCR=8314.36D0) * Molecular weight of dry air PARAMETER (DMD=28.966D0) * Molecular weight of water vapour PARAMETER (DMW=18.016D0) * Mean Earth radius (metre) PARAMETER (S=6378120D0) * Exponent of temperature dependence of water vapour pressure PARAMETER (GAMMA=18.36D0) * Height of tropopause (metre) PARAMETER (HT=11000D0) * Upper limit for refractive effects (metre) PARAMETER (HS=80000D0) INTEGER IN,IS,K,ISTART,I,J DOUBLE PRECISION WLSQ,GB,Z1,A(10),PW0,R0,T0O,DN0,DNDR0, : SK0,F0,RT,TT,DNT,DNDRT,ZT,FT,DNTS,DNDRTS, : SINE,ZTS,FTS,RS,DNS,DNDRS,ZS,FS,REFO, : FE,FO,H,FB,FF,STEP,Z,R,RG,TG, : T,DN,DNDR,F,REFP,REFT * The refraction integrand DOUBLE PRECISION REFI REFI(R,DN,DNDR) = R*DNDR/(DN+R*DNDR) * Set up model atmosphere parameters defined at the observer WLSQ = WL*WL GB= 9.784D0*(1D0-0.0026D0*COS(2D0*PH)-0.00000028D0*H0) Z1 = ((287.604D0+1.6288D0/WLSQ+0.0136D0/(WLSQ*WLSQ)) : *273.15D0/1013.25D0)*1D-6 A(1) = ABS(AS) A(2) = (GB*DMD)/GCR A(3) = A(2)/A(1) A(4) = GAMMA PW0 = UPS*(T0/247.1D0)**A(4) A(5) = PW0*(1D0-DMW/DMD)*A(3)/(A(4)-A(3)) A(6) = P0+A(5) A(7) = Z1*A(6)/T0 A(8) = (Z1*A(5)+11.2684D-6*PW0)/T0 A(9) = (A(3)-1D0)*A(1)*A(7)/T0 A(10) = (A(4)-1D0)*A(1)*A(8)/T0 A(3) = A(3)-2D0 A(4) = A(4)-2D0 * At the observer R0 = S+H0 CALL sla_$ATMT(R0,T0,A,R0,T0O,DN0,DNDR0) SK0 = DN0*R0*SIN(Z0) F0 = REFI(R0,DN0,DNDR0) * At the tropopause in the troposphere RT = S+HT CALL sla_$ATMT(R0,T0,A,RT,TT,DNT,DNDRT) SINE = SK0/(RT*DNT) ZT = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0))) FT = REFI(RT,DNT,DNDRT) * At the tropopause in the stratosphere CALL sla_$ATMS(RT,TT,DNT,A(2),RT,DNTS,DNDRTS) SINE = SK0/(RT*DNTS) ZTS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0))) FTS = REFI(RT,DNTS,DNDRTS) * At the stratosphere limit RS = S+HS CALL sla_$ATMS(RT,TT,DNT,A(2),RS,DNS,DNDRS) SINE = SK0/(RS*DNS) ZS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0))) FS = REFI(RS,DNS,DNDRS) * * Integrate the refraction integral in two parts; first in the * troposphere (K=1), then in the stratosphere (K=2) * REFO = -999.999D0 IS = 16 DO K = 1,2 ISTART = 0 FE = 0D0 FO = 0D0 IF (K.EQ.1) THEN H = (ZT-Z0)/DBLE(IS) FB = F0 FF = FT ELSE H = (ZS-ZTS)/DBLE(IS) FB = FTS FF = FS ENDIF IN = IS-1 IS = IS/2 STEP = H * * Start of iteration loop (terminates when specified precision reached) 200 CONTINUE DO I=1,IN IF (I.EQ.1.AND.K.EQ.1) THEN Z = Z0 + H R = R0 ELSE IF (I.EQ.1.AND.K.EQ.2) THEN Z = ZTS + H R = RT ELSE Z = Z + STEP ENDIF * Given the zenith distance (Z) find R RG = R DO J=1,4 IF (K.EQ.1) THEN CALL sla_$ATMT(R0,T0,A,RG,TG,DN,DNDR) ELSE CALL sla_$ATMS(RT,TT,DNT,A(2),RG,DN,DNDR) ENDIF IF (Z.GT.1D-20) : RG = RG-((RG*DN-SK0/SIN(Z))/(DN+RG*DNDR)) ENDDO R = RG * Find refractive index and integrand at R IF (K.EQ.1) THEN CALL sla_$ATMT(R0,T0,A,R,T,DN,DNDR) ELSE CALL sla_$ATMS(RT,TT,DNT,A(2),R,DN,DNDR) ENDIF F = REFI(R,DN,DNDR) IF (ISTART.EQ.0.AND.MOD(I,2).EQ.0) THEN FE = FE+F ELSE FO = FO+F ENDIF ENDDO * Evaluate the integrand using Simpson's Rule REFP = H*(FB+4D0*FO+2D0*FE+FF)/3D0 * Skip if required precision reached IF (ABS(REFP-REFO).LE.0.5D0*EPS) GO TO 900 * Prepare for next iteration IS = 2*IS IN = IS STEP = H H = H/2D0 FE = FE+FO FO = 0D0 REFO = REFP IF (ISTART.EQ.0) ISTART=1 * Iterate GO TO 200 * Save troposphere component 900 CONTINUE IF (K.EQ.1) REFT=REFP ENDDO * Result REF=REFT+REFP END SUBROUTINE sla_$ATMT (R0,T0,A,R,T,DN,DNDR) *+ * * - - - - - * A T M T * - - - - - * * * Internal routine used by REFRO * * Model atmosphere for the troposphere * * * Given: * RT dp height of tropopause from centre of the Earth (metre) * R0 dp height of observer from centre of the Earth (metre) * T0 dp temperature at the observer (deg K) * A dp(10) constants defined at the observer * R dp current distance from the centre of the Earth (metre) * * Returned: * T dp temperature at R (deg K) * DN dp refractive index at R * DNDR dp rate the refractive index is changing at R * * This routine is a version of the ATMOSTRO routine (C.Hohenkerk, * HMNAO), with trivial modifications. * * Original version by C.Hohenkerk, HMNAO, August 1984 * Starlink version by P.T.Wallace, July 1986 * *+ IMPLICIT NONE DOUBLE PRECISION R0,T0,A(10),R,T,DN,DNDR DOUBLE PRECISION TT0,TT01,TT02 T = T0-A(1)*(R-R0) TT0 = T/T0 TT01 = TT0**(A(3)) TT02 = TT0**(A(4)) DN = 1D0+(A(7)*TT01-A(8)*TT02)*TT0 DNDR = -A(9)*TT01+A(10)*TT02 END SUBROUTINE sla_$ATMS (RT,TT,DNT,GMDR,R,DN,DNDR) *+ * * - - - - - * A T M S * - - - - - * * * Internal routine used by REFRO * * Model atmosphere for the stratosphere * * * Given: * RT dp height of tropopause from centre of the Earth (metre) * TT dp temperature at the tropopause (deg K) * DNT dp refractive index at the tropopause * GMDR dp constant of the atmospheric model = G*MD/R * R dp current distance from the centre of the Earth (metre) * * Returned: * DN dp refractive index at R * DNDR dp rate the refractive index is changing at R * * This routine is a version of the ATMOSSTR routine (C.Hohenkerk, * HMNAO), with trivial modifications. * * Original version by C.Hohenkerk, HMNAO, August 1984 * Starlink version by P.T.Wallace, July 1986 * *+ IMPLICIT NONE DOUBLE PRECISION RT,TT,DNT,GMDR,R,DN,DNDR DOUBLE PRECISION B B = GMDR/TT DN = 1D0+(DNT-1D0)*EXP(-B*(R-RT)) DNDR = -B*(DNT-1D0)*EXP(-B*(R-RT)) END SUBROUTINE sla_REFZ (ZU,REFA,REFB,ZR) *+ * * - - - - - * R E F Z * - - - - - * * * Adjust an unrefracted zenith distance to include the effect of * atmospheric refraction, using the simple A tan Z + B tan**3 Z * model. * * * Given: * ZU dp unrefracted zenith distance of the source (radian) * REFA dp tan Z coefficient (radian) * REFB dp tan**3 Z coefficient (radian) * * Returned: * ZR dp refracted zenith distance (radian) * * Note that this routine applies the adjustment for refraction in * the opposite sense to the usual one - it takes an unrefracted * (in vacuo) position and produces an observed (refracted) * position, whereas the basic A tan Z + B tan**3 Z model strictly * applies to the case where a refracted position is available and * must be corrected for refraction. This requires an inverted form of * the refraction expression; the formula used here is taken from * RGO/LaPalma Technical Note No. 42, "A control system for * the William Herschel Telescope", by R.Laing. * * For ZD > 85 deg, a fixed value ZD = 85 deg is used. * * See also the routine sla_REFV, which performs the adjustment in * Cartesian Az/El coordinates. * * P.T.Wallace Starlink July 1986 * *+ IMPLICIT NONE DOUBLE PRECISION ZU,REFA,REFB,ZR DOUBLE PRECISION ZL,S,C,T,TSQ,TCU * Keep ZD within usable range ZL = MIN(ABS(ZU),1.483529864D0) * Functions of ZD S = SIN(ZL) C = COS(ZL) T = S/C TSQ = T*T TCU = T*TSQ * Refracted ZD ZR = ZU-(REFA*T+REFB*TCU)/(1D0+(REFA+3D0*REFB*TSQ)/(C*C)) END SUBROUTINE sla_REFCO (H0,T0,P0,UPS,WL,PH,AS,EPS,REFA,REFB) *+ * * - - - - - - * R E F C O * - - - - - - * * * Determine constants A and B in atmospheric refraction model * dZ = A tan Z + B tan**3 Z (for optical wavelengths). * * * Given: * H0 dp height of the observer above sea level (metre) * T0 dp ambient temperature at the observer (deg K) * P0 dp pressure at the observer (millibar) * UPS dp relative humidity at the observer (range 0-1) * WL dp effective wavelength of the source (micrometre) * PH dp latitude of the observer (radian, astronomical) * AS dp temperature lapse rate in the troposphere (degK/metre) * EPS dp precision required to terminate iteration (radian) * * Returned: * REFA dp tan Z coefficient (radian) * REFB dp tan**3 coefficient (radian) * * Called: sla_REFRO * * Typical values for the AS and EPS arguments would be 0.0065D0 and * 1D-10 respectively. * * Relative to the comprehensive refraction model used by this routine, * the simple A tan Z + B tan**3 Z formula achieves 0.5 arcsec accuracy * for ZD < 80 deg, 0.01 arcsec accuracy for ZD < 60 deg, and * 0.001 arcsec accuracy for ZD < 45 deg. * * P.T.Wallace Starlink July 1987 * *+ IMPLICIT NONE DOUBLE PRECISION H0,T0,P0,UPS,WL,PH,AS,EPS,REFA,REFB DOUBLE PRECISION ATN1,ATN4,R1,R2 * Sample zenith distances: arctan(1) and arctan(4) PARAMETER (ATN1=0.7853981634D0, : ATN4=1.325817664D0) * Determine refraction for the two sample zenith distances CALL sla_REFRO(ATN1,H0,T0,P0,UPS,WL,PH,AS,EPS,R1) CALL sla_REFRO(ATN4,H0,T0,P0,UPS,WL,PH,AS,EPS,R2) * Solve for refraction constants REFA = (64D0*R1-R2)/60D0 REFB = (R2-4D0*R1)/60D0 END SUBROUTINE RDCAL(X,Y,PCRA,PCDEC,ITEL,RA,DEC) * ADJUST IDEAL PLATE COORDINATES (X,Y) ACCORDING TO * TELESCOPE TYPE GIVING TANGENTIAL COORDS (XI,ETA) * THEN CALCULATES (RA,DEC) * GIVEN: * X IDEAL PLATE COORDINATE * Y IDEAL PLATE COORDINATE * PCRA PLATE CENTRE RA * PCDEC PLATE CENTRE DEC * ITEL TELESCOPE TYPE * RETURNED: * RA CALCULATED RA * DEC CALCULATED DEC IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION X,Y,PCRA,PCDEC INTEGER ITEL DOUBLE PRECISION RA,DEC * ADJUST FOR TELESCOPE TYPE GO TO (11,12,13,14),ITEL * ASTROGRAPH 11 CONTINUE XI=X ETA=Y GO TO 20 * SCHMIDT 12 CONTINUE CALL SCHMIT(X,Y,CURV) XI=X*CURV ETA=Y*CURV GO TO 20 * AAT DOUBLET 13 CONTINUE XI=X ETA=Y CALL PNCUSH(ITEL,XI,ETA,'DU') GO TO 20 * AAT TRIPLET 14 CONTINUE XI=X ETA=Y CALL PNCUSH(ITEL,XI,ETA,'DU') 20 CONTINUE * COMPUTE (RA,DEC) CALL DTP2S(XI,ETA,PCRA,PCDEC,RA,DEC) END SUBROUTINE PNCUSH(ITEL,X,Y,KDIRN) * ALLOW FOR PINCUSHION DISTORTION * GIVEN: * ITEL TELESCOPE TYPE : 3 = AAT2, 4 = AAT3 * X,Y PLATE COORDINATES (RADIANS) * KDIRN DIRECTION OF ADJUSTMENT * RETURNED: * X,Y ADJUSTED * KDIRN VALUES: * 'UD' FROM TANGENT-PLANE COORDINATES * PREDICTS DISTORTED COORDINATES * 'DU' FROM DISTORTED COORDINATES * PREDICTS TANGENT-PLANE COORDINATES IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER ITEL DOUBLE PRECISION X,Y CHARACTER*2 KDIRN * PINCUSHION DISTORTION COEFFICIENTS : AAT DOUBLET, TRIPLET * (DOUBLET USED AS ERROR DEFAULT) PARAMETER (PCAAT2=147.069D0) PARAMETER (PCAAT3=178.585D0) * PICK UP CORRECT DISTORTION COEFFICIENT D=PCAAT3 IF(ITEL.NE.4) D=PCAAT2 * DISTANCE FROM CENTRE RSQ=X*X + Y*Y * USEFUL EXPRESSION DRSQ=D*RSQ * * DEFAULT TO NO DISTORTION IF ILLEGAL KDIRN RFACTR=1D0 * UNDISTORTED TO DISTORTED IF(KDIRN.EQ.'UD') RFACTR=1D0+DRSQ * DISTORTED TO UNDISTORTED IF(KDIRN.EQ.'DU') RFACTR=(2D0*DRSQ+1D0)/(3D0*DRSQ+1D0) * ADJUST X, Y X=X*RFACTR Y=Y*RFACTR END SUBROUTINE SCHMIT(A,B,CURV) * CALCULATE CURVATURE CORRECTION FOR SCHMIDT PLATES * X = XI / CURV , Y = ETA / CURV * XI = X * CURV , ETA = Y * CURV * GIVEN: * A TANGENTIAL COORD (XI) OR IDEAL PLATE COORD (X) * B TANGENTIAL COORD (ETA) OR IDEAL PLATE COORD (Y) * RETURNED: * CURV CURVATURE CORRECTION FACTOR IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A,B,CURV R=SQRT(A*A + B*B) CURV=1D0 IF(ABS(R).GT.1D-20) CURV=(SIN(R)/COS(R))/R END SUBROUTINE RMS(LUR,LUS,SIG,NREF) * CALCULATE AND PRINTS RMS OF RESIDUALS * DR, DD, EW, SKY * GIVEN: * LUR LOGICAL UNIT FOR FULL REPORT * LUS LOGICAL UNIT FOR SYNOPSIS * NREF NUMBER OF REFERENCE STARS * MODIFIED: * SIG RMS SUMMATION VECTOR IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER LUR,LUS,NREF DOUBLE PRECISION SIG(4) * ARC SECONDS TO RADIANS PARAMETER (AS2R=4.8481368110953599D-6) REFN=DBLE(NREF) DO I=1,4 SIG(I)=SQRT(SIG(I)/REFN)/AS2R END DO WRITE(LUR,1000) SIG 1000 FORMAT(/75X,'RMS :',F9.2,2F10.2,F9.2/) WRITE(LUS,1001) SIG(3),SIG(2),SIG(4) 1001 FORMAT(/17X,'RMS :',F9.2,2F10.2) END SUBROUTINE DFLTIN(STRING,NSTRT,LENFLD,LENSTR,DRESLT,JFLAG) * * * CONVERT FREE-FORMAT INPUT INTO DOUBLE PRECISION FLOATING POINT * * GIVEN: * STRING CHARACTER STRING CONTAINING FIELD TO BE DECODED * NSTRT POINTER TO 1ST CHARACTER OF FIELD IN STRING * LENFLD FIELD LENGTH (CHARACTERS) * LENSTR STRING LENGTH (CHARACTERS) * * RETURNED: * NSTRT INCREMENTED * DRESLT DOUBLE PRECISION RESULT * JFLAG -1 = -OK, 0 = +OK, 1 = NULL FIELD, 2 = ERROR * * * NOTES: * 1 THE BASIC FORMAT IS #^.^@#^ WHERE # MEANS + OR -, * ^ MEANS A DECIMAL SUBFIELD AND @ MEANS D OR E. * 2 BLANKS: * LEADING BLANKS ARE IGNORED. * EMBEDDED BLANKS ARE ALLOWED ONLY AFTER # AND D OR E, AND * AFTER . WHERE THE FIRST ^ IS ABSENT. * TRAILING BLANKS ARE IGNORED; THE FIRST SIGNIFIES * END OF DECODING AND SUBSEQUENT ONES ARE SKIPPED. * 3 FIELD SEPARATORS: * ANY CHARACTER OTHER THAN +,-,0-9,.,D,E OR BLANK MAY BE * USED TO END A FIELD. COMMA IS RECOGNISED BY DFLTIN * AS A SPECIAL CASE; IT IS SKIPPED, LEAVING THE * POINTER ON THE NEXT CHARACTER. SEE 11, BELOW. * 4 BOTH SIGNS ARE OPTIONAL. THE DEFAULT IS +. * 5 THE MANTISSA DEFAULTS TO 1. * 6 THE EXPONENT DEFAULTS TO D0. * 7 THE DECIMAL SUBFIELDS MAY BE OF ANY LENGTH. * 8 THE DECIMAL POINT IS OPTIONAL FOR INTEGERS. * 9 A NULL FIELD IS ONE THAT DOES NOT BEGIN WITH * +,-,0-9,.,D OR E, OR CONSISTS ENTIRELY OF BLANKS. * IF THE FIELD IS NULL, JFLAG IS SET TO 1 AND * RESULT IS LEFT UNTOUCHED. * 10 NSTRT = 1 FOR THE FIRST CHARACTER IN THE STRING. * 11 ON RETURN FROM DFLTIN, NSTRT POINTS TO THE CHARACTER * FOLLOWING THE FINAL CHARACTER OF THE FIELD. IF A * SEPARATOR IS BEING USED - OTHER THAN COMMA - NSTRT * MUST BE INCREMENTED BEFORE THE NEXT CALL TO DFLTIN. * 12 ERRORS (JFLAG=2) OCCUR WHEN: * A) A +, -, D OR E IS LEFT UNSATISFIED. * B) THE DECIMAL POINT IS PRESENT WITHOUT AT LEAST * ONE DECIMAL SUBFIELD. * C) AN EXPONENT MORE THAN 100 HAS BEEN PRESENTED. * VARIOUS ARITHMETIC FAULTS ARE POSSIBLE WHICH WILL * BE REPORTED BY THE ENVIRONMENT * RATHER THAN BY DFLTIN. * 13 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 * UNTIL THE WHOLE FIELD HAS BEEN READ. * 14 CERTAIN HIGHLY UNLIKELY COMBINATIONS OF MANTISSA & * EXPONENT CAN CAUSE ARITHMETIC FAULTS DURING THE * DECODE, DESPITE THE FACT THAT THEY TOGETHER * COULD BE CONSTRUED AS A VALID NUMBER. * 15 DECODING IS LEFT TO RIGHT, ONE PASS. * 16 END OF FIELD MAY OCCUR IN ONE OF THREE WAYS: * A) AS DICTATED BY THE FIELD LENGTH, LENFLD. * B) AS DICTATED BY THE STRING LENGTH, LENSTR. * C) DETECTED DURING THE DECODE. * (C OVERRIDES B, AND B OVERRIDES A.) IMPLICIT DOUBLE PRECISION (D) CHARACTER*(*) STRING INTEGER NSTRT,LENFLD,LENSTR DOUBLE PRECISION DRESLT INTEGER JFLAG DATA DTEN,DTEN10/1D1,1D10/ * SET UP POINTERS: * CURRENT CHARACTER JPTR=NSTRT * END OF FIELD (OR STRING) JNXFLD=MIN0(JPTR+LENSTR-1,LENFLD) * 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 IDCHF(STRING,JPTR,JNXFLD,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 IDCHF(STRING,JPTR,JNXFLD,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*DTEN+DIGIT CALL IDCHF(STRING,JPTR,JNXFLD,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 IDCHF(STRING,JPTR,JNXFLD,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 IDCHF(STRING,JPTR,JNXFLD,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*DTEN+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 IDCHF(STRING,JPTR,JNXFLD,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 IDCHF(STRING,JPTR,JNXFLD,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 IDCHF(STRING,JPTR,JNXFLD,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 JPTR=JPTR-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*DTEN10 NEXP=NEXP-10 GO TO 1400 1410 CONTINUE IF(NEXP.LT.1) GO TO 1600 DMANT=DMANT*DTEN NEXP=NEXP-1 GO TO 1410 * NEGATIVE EXPONENT: SCALE DOWN 1500 CONTINUE IF(NEXP.GT.-10) GO TO 1510 DMANT=DMANT/DTEN10 NEXP=NEXP+10 GO TO 1500 1510 CONTINUE IF(NEXP.GT.-1) GO TO 1600 DMANT=DMANT/DTEN 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 IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT) GO TO (1720, 1620, 1720, 1720, 1720, 1720, 9900, 1720, 9900),NVEC * 0-9 SP D/E . + - , 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 IDCHF(STRING,JPTR,JNXFLD,NVEC,NDIGIT,DIGIT) * * IDENTIFY NEXT CHARACTER IN STRING * * GIVEN: * STRING CHARACTER STRING * JPTR POINTER TO CHARACTER TO BE IDENTIFIED * JNXFLD POINTER TO END OF FIELD * * RETURNED: * JPTR INCREMENTED UNLESS END OF FILED * NVEC VECTOR FOR IDENTIFIED CHARACTER * NDIGIT INTEGER 0-9 IF CHARACTER WAS A NUMERAL * DIGIT DOUBLE PRECISION EQUIVALENT OF NDIGIT * * NVEC TAKES THE FOLLOWING VALUES: * 1 0-9 * 2 SPACE * 3 D,E * 4 . * 5 + * 6 - * 7 , * 8 ELSE * 9 END OF FIELD CHARACTER*(*) STRING INTEGER JPTR,JNXFLD,NVEC,NDIGIT DOUBLE PRECISION DIGIT * CHARACTER/VECTOR TABLES CHARACTER KCTAB(17) INTEGER KVTAB(17) DATA KCTAB/'0','1','2','3','4','5','6','7','8','9', * ' ','D','E','.','+','-',','/ DATA KVTAB/10*1,2,2*3,4,5,6,7/ * HANDLE END OF FIELD IF(JPTR.GT.JNXFLD) THEN NVEC=9 ELSE * NOT END OF FIELD * IDENTIFY CHARACTER DO NCHAR=1,17 IF(KCTAB(NCHAR).EQ.STRING(JPTR:JPTR)) GO TO 2200 END DO * UNRECOGNISED NVEC=8 GO TO 2300 * RECOGNISED 2200 CONTINUE NVEC=KVTAB(NCHAR) * ALLOW FOR NUMERALS NDIGIT=NCHAR-1 DIGIT=DBLE(FLOAT(NDIGIT)) * INCREMENT POINTER 2300 CONTINUE JPTR=JPTR+1 END IF * RETURN END SUBROUTINE DMATRX(N,A,Y,D,JF,W,IW) * * MATRIX INVERSION & SOLUTION TO SIMULTANEOUS EQUATIONS * * DOUBLE PRECISION * * FOR THE SET OF N SIMULTANEOUS EQUATIONS IN N UNKNOWNS: * A.Y = X * WHERE: * A IS A NON-SINGULAR N X N MATRIX * Y IS THE VECTOR OF N UNKNOWNS * X IS THE KNOWN VECTOR * DMATRX COMPUTES: * THE INVERSE OF MATRIX A * THE DETERMINANT OF MATRIX A * THE VECTOR OF N UNKNOWNS * * ARGUMENTS: * SYMBOL TYPE DIMENSION BEFORE AFTER * N I NO. OF UNKNOWNS UNCHANGED * A DP (N,N) MATRIX INVERSE * Y DP (N) VECTOR SOLUTION * D DP - DETERMINANT * * JF I - SINGULARITY FLAG * W DP (N,N) - WORKSPACE * IW I (N) - WORKSPACE * * * JF IS THE SINGULARITY FLAG. IF THE MATRIX IS NON-SINGULAR, * JF=0 IS RETURNED. IF THE MATRIX IS SINGULAR, JF=-1 & D=0.0 ARE * RETURNED. IN THE LATTER CASE, THE CONTENTS OF ARRAY A ON RETURN * ARE UNDEFINED. * * PEDIGREE: * BY PAUL MURDIN OUT OF UNIVERSITY OF ROCHESTER, NEW YORK. * * ALGORITHM: * NO DESCRIPTION AVAILABLE. SEEMS TO BE BASED ON A STANDARD * METHOD BUT WITH RIGOUR IN THE SELECTION OF PIVOT POINTS * COMPROMISED TO IMPROVE SPEED AND SIZE AT THE EXPENSE OF ACCURACY. * * SPEED: * VERY FAST * * ACCURACY: * FAIRLY ACCURATE - ERRORS 1 TO 4 TIMES THOSE OF ROUTINES OPTIMISED * FOR ACCURACY. * DOUBLE PRECISION A(N,N),Y(N),D,W(N,N) INTEGER IW(N) DOUBLE PRECISION AMX,T,AKK,YK,AIK,SFA PARAMETER (SFA=1D-20) JF=0 D=1D0 DO K=1,N AMX=DABS(A(K,K)) IMX=K IF (K.NE.N) THEN DO I=K+1,N T=DABS(A(I,K)) IF (T.GT.AMX) THEN AMX=T IMX=I END IF END DO END IF IF (AMX.LT.SFA) THEN JF=-1 ELSE IF (IMX.NE.K) THEN DO J=1,N T=A(K,J) A(K,J)=A(IMX,J) A(IMX,J)=T END DO T=Y(K) Y(K)=Y(IMX) Y(IMX)=T D=-D END IF IW(K)=IMX AKK=A(K,K) D=D*AKK IF (DABS(D).LT.SFA) THEN JF=-1 ELSE AKK=1D0/AKK A(K,K)=AKK DO J=1,N IF (J.NE.K) A(K,J)=A(K,J)*AKK END DO YK=Y(K)*AKK Y(K)=YK DO I=1,N AIK=A(I,K) IF (I.NE.K) THEN DO J=1,N IF (J.NE.K) A(I,J)=A(I,J)-AIK*A(K,J) END DO Y(I)=Y(I)-AIK*YK END IF END DO DO I=1,N IF (I.NE.K) A(I,K)=-A(I,K)*AKK END DO END IF END IF END DO IF (JF.NE.0) THEN D=0D0 ELSE DO K=1,N NP1MK=N+1-K KI=IW(NP1MK) IF (NP1MK.NE.KI) THEN DO I=1,N T=A(I,NP1MK) A(I,NP1MK)=A(I,KI) A(I,KI)=T END DO END IF END DO END IF END SUBROUTINE DR2AF(NDP,ANGLE,SIGN,IDMSF) * * * CONVERT ANGLE FOR OUTPUT * ------------------------ * * * GIVEN: * NDP NUMBER OF DECIMAL PLACES OF ARC SECONDS * ANGLE ANGLE IN RADIANS (DOUBLE PRECISION) * * RETURNED: * SIGN '+' OR '-' * IDMSF DEGREES, ARC MINUTES, ARC SECONDS, FRACTION * * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * INTEGER NDP DOUBLE PRECISION ANGLE CHARACTER*1 SIGN INTEGER IDMSF(4) * RADIANS TO ARC SECONDS DOUBLE PRECISION R2AS PARAMETER (R2AS=206264.80624709636D0) * HANDLE SIGN SIGN='+' IF(ANGLE.LT.0D0) SIGN='-' * FIELD UNITS IN TERMS OF LEAST SIGNIFICANT FIGURE RAS=REAL(10**MAX(NDP,0)) RAM=RAS*60D0 RD=RAM*60D0 * ROUND ANGLE AND EXPRESS IN SMALLEST UNITS REQUIRED A=NINT(RAS*R2AS*ABS(ANGLE)) * SEPARATE INTO FIELDS AD=INT(A/RD) A=A-AD*RD AM=INT(A/RAM) A=A-AM*RAM AS=INT(A/RAS) AF=A-AS*RAS * RETURN RESULTS IDMSF(1)=INT(AD) IDMSF(2)=INT(AM) IDMSF(3)=INT(AS) IDMSF(4)=INT(AF) END SUBROUTINE DR2TF(NDP,ANGLE,SIGN,IHMSF) * * * CONVERT ANGLE FOR OUTPUT, EXPRESSED AS TIME * ------------------------------------------- * * * GIVEN: * NDP NUMBER OF DECIMAL PLACES OF SECONDS * ANGLE ANGLE IN RADIANS (DOUBLE PRECISION) * * RETURNED: * SIGN '+' OR '-' * IHMSF HOURS, MINUTES, SECONDS, FRACTION * IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NDP DOUBLE PRECISION ANGLE CHARACTER*1 SIGN INTEGER IHMSF(4) * RADIANS TO SECONDS DOUBLE PRECISION R2S PARAMETER (R2S=13750.987083139757D0) * HANDLE SIGN SIGN='+' IF(ANGLE.LT.0D0) SIGN='-' * FIELD UNITS IN TERMS OF LEAST SIGNIFICANT FIGURE RS=REAL(10**MAX(NDP,0)) RM=RS*60D0 RH=RM*60D0 * ROUND ANGLE AND EXPRESS IN SMALLEST UNITS REQUIRED A=NINT(RS*R2S*ABS(ANGLE)) * SEPARATE INTO FIELDS AH=INT(A/RH) A=A-AH*RH AM=INT(A/RM) A=A-AM*RM AS=INT(A/RS) AF=A-AS*RS * RETURN RESULTS IHMSF(1)=INT(AH) IHMSF(2)=INT(AM) IHMSF(3)=INT(AS) IHMSF(4)=INT(AF) END DOUBLE PRECISION FUNCTION DRANGE(ANGLE) * * * NORMALISE ANGLE INTO RANGE +/- PI * --------------------------------- * * * ANGLE IS THE ANGLE IN RADIANS * BOTH ANGLE AND THE RESULT ARE DOUBLE PRECISION * DOUBLE PRECISION ANGLE DOUBLE PRECISION DPI,D2PI PARAMETER (DPI=3.141592653589793238462643D0) PARAMETER (D2PI=6.283185307179586476925287D0) DRANGE=DMOD(ANGLE,D2PI) IF(DABS(DRANGE).GE.DPI) DRANGE=DRANGE-DSIGN(D2PI,ANGLE) END DOUBLE PRECISION FUNCTION DRANRM(ANGLE) * * * NORMALISE ANGLE INTO RANGE 0-2PI * -------------------------------- * * * ANGLE IS THE ANGLE IN RADIANS * BOTH ANGLE AND THE RESULT ARE DOUBLE PRECISION DOUBLE PRECISION ANGLE DOUBLE PRECISION D2PI PARAMETER (D2PI=6.283185307179586476925287D0) DRANRM=DMOD(ANGLE,D2PI) IF(DRANRM.LT.0D0) DRANRM=DRANRM+D2PI END SUBROUTINE DS2TP(RA,DEC,RAZ,DECZ,XI,ETA) * * * PROJECTION OF SPHERICAL COORDINATES ONTO TANGENT PLANE * ('GNOMONIC' PROJECTION - 'STANDARD COORDINATES') * * * GIVEN: * RA,DEC SPHERICAL COORDINATES OF POINT TO BE PROJECTED * RAZ,DECZ SPHERICAL COORDINATES OF TANGENT POINT * * RETURNED: * XI,ETA RECTANGULAR COORDINATES ON TANGENT PLANE * IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION RA,DEC,RAZ,DECZ,XI,ETA SDECZ=SIN(DECZ) SDEC=SIN(DEC) CDECZ=COS(DECZ) CDEC=COS(DEC) RADIF=RA-RAZ SRADIF=SIN(RADIF) CRADIF=COS(RADIF) DENOM=SDEC*SDECZ+CDEC*CDECZ*CRADIF XI=CDEC*SRADIF/DENOM ETA=(SDEC*CDECZ-CDEC*SDECZ*CRADIF)/DENOM END SUBROUTINE DTP2S(XI,ETA,RAZ,DECZ,RA,DEC) * * * TRANSFORM TANGENT PLANE COORDINATES INTO SPHERICAL * * * GIVEN: * XI,ETA TANGENT PLANE RECTANGULAR COORDINATES * RAZ,DECZ SPHERICAL COORDINATES OF TANGENT POINT * * RETURNED: * RA,DEC SPHERICAL COORDINATES * IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION XI,ETA,RAZ,DECZ,RA,DEC SDECZ=SIN(DECZ) CDECZ=COS(DECZ) DENOM=CDECZ-ETA*SDECZ RADIF=ATAN2(XI,DENOM) RA=DRANRM(RADIF+RAZ) DEC=ATAN2(SDECZ+ETA*CDECZ,SQRT(XI*XI+DENOM*DENOM)) END SUBROUTINE MMP(RM0, DM0, PR, PD, EPOCH0, EPOCH1, RM1, DM1) * * GIVEN: * RM0, DM0 MEAN RA, DEC AT 1ST EPOCH * PR, PD PROPER MOTIONS IN RA, DEC AT SPECIFIED EPOCH * (RADIANS PER YEAR) * EPOCH0 1ST EPOCH * EPOCH1 2ND EPOCH * * RETURNED: * RM1, DM1 MEAN RA, DEC AT 2ND EPOCH * IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION RM0, DM0, PR, PD, EPOCH0, EPOCH1, RM1, DM1 DOUBLE PRECISION RMATP(3,3), V(3) * PROPER MOTIONS TAU= EPOCH1 - EPOCH0 R= RM0 + PR*TAU D= DM0 + PD*TAU * PRECESSION MATRIX CALL PREC(EPOCH0, EPOCH1, RMATP) * MEAN TO MEAN CALL DCS2C(R, D, V) CALL RMATMU(V, RMATP) CALL DCC2S(R, DM1, V) RM1= DRANRM(R) END SUBROUTINE PREC(EPOCH0, EPOCH1, S) * * GIVEN: * EPOCH0 STARTING EPOCH (YEARS AD) * EPOCH1 ENDING EPOCH (YEARS AD) * * RETURNED: * S PRECESSION MATRIX * IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION EPOCH0, EPOCH1, S(3,3) TSER(C0,C1,C2)=DRANGE(C0+(C1+C2*T)*T) * TIME FROM 1900.0 TO ORIGINAL EPOCH (TROPICAL CENTURIES) T0=(EPOCH0 - 1900D0)/100D0 * TIME SINCE ORIGINAL EPOCH (TROPICAL CENTURIES) T= (EPOCH1 - EPOCH0)/100D0 TSQ=T*T * PRECESSION PARAMETERS Z0=TSER(0., 6.768D-6*T0 + 1.11713192D-2, 1.464D-6) Z=Z0 + 3.835D-6*TSQ TH=TSER(0D0, -4.135D-6*T0 + 9.7189726D-3, -2.065D-6) SZ0=SIN(Z0) CZ0=COS(Z0) SZ=SIN(Z) CZ=COS(Z) STH=SIN(TH) CTH=COS(TH) CTHSZ=CTH*SZ CTHCZ=CTH*CZ * ROTATION MATRIX S(1,1)= CZ0*CTHCZ - SZ0*SZ S(1,2)=-SZ0*CTHCZ - CZ0*SZ S(1,3)=-STH*CZ S(2,1)= CZ0*CTHSZ + SZ0*CZ S(2,2)=-SZ0*CTHSZ + CZ0*CZ S(2,3)=-STH*SZ S(3,1)= CZ0*STH S(3,2)=-SZ0*STH S(3,3)= CTH END SUBROUTINE PMC(PR, PD, EPOCH, DATE, DR, DD) * * GIVEN: * PR, PD PROPER MOTIONS IN RA, DEC AT SPECIFIED EPOCH * (RADIANS PER YEAR) * EPOCH YEAR AD OF MEAN EQUATOR & EQUINOX * DATE DAYS SINCE JD 2400000.0 FOR APPARENT PLACE * * RETURNED: * DR, DD TOTAL PROPER MOTION CORRECTIONS (RADIANS) * IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION PR, PD, EPOCH, DATE, DR, DD * TIME SINCE ORIGINAL EPOCH (TROPICAL YEARS) TAU= (DATE - 15020.313D0)/365.2421957D0 - (EPOCH - 1900D0) * PROPER MOTION CORRECTIONS (RA,DEC VECTOR) DR= PR*TAU DD= PD*TAU END SUBROUTINE DCS2C(A,B,V) * * A IS RA/HA/LONG/AZ * B IS DEC/LAT/ZD * V IS X,Y,Z * IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A, B, V(3) AR= DRANGE(A) BR= DRANGE(B) COSBR= COS(BR) V(1)= COS(AR)*COSBR V(2)= SIN(AR)*COSBR V(3)= SIN(BR) END SUBROUTINE DCC2S(A,B,V) * * A IS RA/HA/LONG/AZ * B IS DEC/LAT/ZD * V IS X,Y,Z * IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A, B, V(3) X=V(1) Y=V(2) A=ATAN2(Y,X) B=ATAN2(V(3), SQRT(X*X + Y*Y)) END SUBROUTINE RMATMU(A,B) * * PERFORMS THE 3-D FORWARD UNITARY TRANSFORMATION: * VECTOR A = MATRIX B * VECTOR A * DOUBLE PRECISION A(3), B(3,3) DOUBLE PRECISION C(3) DO J=1,3 C(J)=0D0 DO I=1,3 C(J)=C(J)+A(I)*B(J,I) END DO END DO DO J=1,3 A(J)=C(J) END DO END FUNCTION RJDATE ( NYEAR, NDIY ) * * * JULIAN DATE FROM YEAR & DAY-IN-YEAR * * * ARGUMENTS: * * NYEAR GREGORIAN YEAR AD (INTEGER; MUST BE 1901 THROUGH 2100) * NDIY DAY IN YEAR NUMBER (INTEGER; JAN 1ST = DAY 1) * * RESULT IS THE JULIAN DATE MINUS 2400000.0 FOR 0 HRS UT * RJDATE=REAL(NDIY+(NYEAR-1)/4)+REAL(NYEAR-1901)*365.0+14909.5 END