C ************************************************************************** C Program Findtarget C -------------------------------------------------------------------------- C Date Programmer Version description C ---- ---------- ------------------- C Aug 93 Cox, Lallo V1.0 All modules in one file C C Dec 93 Lallo V1.1 Support of multiple targets & C inclusion of guidestar separation C differences in output file. C C Feb 94 Cox V1.2 Enhanced error checking of input C C June 94 Cox V1.3 Synchronized guide star input and C moved fgs offsets into minischf.dat C file. Restructured minischf.dat. C Produces runtime status messages. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FILE NAME: Findtarget.f C V1.0 C SunOS 4.1.3 C C PURPOSE: C C This program takes OMSAT format files of fgs data (encoder angles, pmt C counts) and net spacecraft velocities and a file of guidestar RA & Dec, C plus a target position in v2v3 or ra and dec, and from this determines C the accurate spacecraft attitude corrected for all current field dis- C tortions in the fgs's, alignment matrices, and velocity abberration. C From the attitude it then maps the supplied target ra and dec onto v2v3 C or vice-versa. C C CALLING SEQUENCE: findtarget C C INPUTS: C C file created by omsat containing the following six hi-res parameters: C 4 pmt counts, 2 encoder angles. These parameters are found under the C FGS/FGE OMSAT option. One file must be created for each fgs. A separate C file may also contain the 3 hi-res net velocity parameters (found C under PCS OUTPUT OMSAT option) and the last input file must contain the C following: Line 1- the words RA or V2 in any location or case, depending C on which coordinate system position you want to know. Line 2- the ra and C dec of one guidestar followed by the word "FGS" and the fgs#. Line 3- C same as line 2 but for the second guidestar. Line 4 though N- the v2v3 or C ra and dec of the targets. C C OUTPUTS: C C file called guidestars.v2v3 containing calculated v2,v3 of the guidestars, C and the file called either target.v2v3 or target.radec depending on which C coordinates were found. C C BUILD INSTRUCTIONS: C C UNIX- %f77 -e -misalign findtarget.f -o findtarget C VMS- $FOR/extend_source findtarget.for C $LINK findtarget C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FOR QUESTIONS CONCERNING THIS CODE, SEE COLIN COX cox@stsci.edu or C MATTHEW LALLO lallo@stsci.edu C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PROGRAM FINDTARGET CALL finelock TYPE *, ' Found guide star V2 V3 values' CALL targetpos END C*************************************************************************** SUBROUTINE FINELOCK C*************************************************************************** C C To derive V2V3 coordinates from FGS data in fine lock C IMPLICIT NONE INTEGER schf REAL*8 alever(3), blever(3), tvs(3,3,3), k1x(3), k1y(3) REAL*8 gm(3), alpha(6,6,3), beta(6,6,3) LOGICAL found, filepresent, first, fine CHARACTER*35 time, time2, result, fes1, fes2, response, guide, mininame CHARACTER*80 firstline1, firstline2 INTEGER*4 fgsA, fgsB, pmtxa1, pmtxb1, pmtya1, pmtyb1 INTEGER*4 pmtxa2, pmtxb2, pmtya2, pmtyb2 REAL*8 rpmtxa1, rpmtxb1, rpmtya1, rpmtyb1 REAL*8 rpmtxa2, rpmtxb2, rpmtya2, rpmtyb2 REAL*8 theta_a1, theta_b1,theta_a2, theta_b2 REAL*8 fgs_ima_vec(3), fgs_obj_vec(3), fgs_pos(3) REAL*8 fgs_veh_vec1(3),fgs_veh_vec2(3), hours, hours2, sec, seconds, seconds2 INTEGER*4 year, day, day0, hour, min, pos, pos1, pos2, pos3 INTEGER*4 stat1, stat2, stat3, dataline1, dataline2 REAL*8 aoffset(3), boffset(3) PARAMETER (schf=4) !Unit number for schf file C DATA aoffset/0.0D0, 0.0D0, 0.0D0/ C DATA aoffset/-187.9348, 2101.9382, -1882.4137/ ! in arcsec Dec 1992 solution C DATA boffset/0.0D0, 0.0D0, 0.0D0/ C Extract necessary data from "MINISCHF" file C Temporary change to allow experimenting with different SCHF files. mininame = 'minischf.dat' INQUIRE(file=mininame,exist=filepresent) If (.not.filepresent) then TYPE *,' THE MINISCHF.DAT FILE WAS NOT FOUND. IT IS USED BY THE' TYPE *,' CODE TO OBTAIN FGS CALIBRATION INFORMATION. A COPY OF' TYPE *,' THIS FILE CAN BE OBTAINED UNDER /tib/u1/lallo/public.' TYPE *,' PLEASE PROVIDE AND THEN RUN FINDTARGET AGAIN.' STOP ELSE OPEN (schf,file= mininame, type='old') ENDIF Call READSCHF(schf, alever, blever,tvs, k1x, k1y, aoffset, boffset, gm, X alpha, beta, found) IF (.NOT. found) STOP CLOSE(schf) TYPE *, ' Is the data in fine lock? ' READ '(A)', response fine = ((response(1:1) .EQ. 'y') .OR. (response(1:1) .EQ. 'Y') ) C Read fine error signal data 27 IF (fine) THEN TYPE *, ' First FGS fine error signal file ' READ '(A)', fes1 INQUIRE(file=fes1,exist=filepresent) IF (filepresent) then OPEN (1, file = fes1,type='old') ELSE type *,' THERE IS NO FILE NAMED ',fes1 type *,' PLEASE RE-ENTER' goto 27 ENDIF 28 TYPE *, ' Second FGS fine error signal file ' READ '(A)', fes2 INQUIRE(file=fes2,exist=filepresent) IF (filepresent) then OPEN (1, file = fes1,type='old') ELSE type *,' THERE IS NO FILE NAMED ',fes2 type *,' PLEASE RE-ENTER' goto 28 ENDIF OPEN (2, file = fes2,type='old') ELSE ! Coarse track 29 TYPE *, ' Guide star file ' READ '(A)', guide INQUIRE(file=guide,exist=filepresent) IF (filepresent) then OPEN (1, file = guide,type='old') ELSE type *,' THERE IS NO FILE NAMED ', guide type *,' PLEASE RE-ENTER' GO TO 29 END IF END IF ! Fine lock result = 'guidestars.v23' OPEN (3, file=result) IF (fine) THEN READ (1, '(//A)'), firstline1 pos=INDEX(firstline1,'F') READ(firstline1,'(X,I1)') fgsA READ (1, '(A)'), firstline1 READ (1, '(A)'), firstline1 READ (2, '(//A)'), firstline2 pos=INDEX(firstline2,'F') READ(firstline2,'(X,I1)') fgsB READ (2, '(A)'), firstline2 READ (2, '(A)'), firstline2 ELSE ! coarse READ (1, '(//A)'), firstline1 pos1 = INDEX(firstline1, 'F1') pos2 = INDEX(firstline1, 'F2') pos3 = INDEX(firstline1, 'F3') IF (pos1 .EQ. 0) THEN IF (pos2 .LT. pos3) THEN fgsA = 2 fgsB = 3 ELSE fgsA = 3 fgsB = 2 END IF END IF IF (pos2 .EQ. 0) THEN IF (pos3 .LT. pos1) THEN fgsA = 3 fgsB = 1 ELSE fgsA = 1 fgsB = 3 END IF END IF IF (pos3 .EQ. 0) THEN IF (pos1 .LT. pos2) THEN fgsA = 1 fgsB = 2 ELSE fgsA = 2 fgsB = 1 END IF END IF TYPE *, ' FGS A =',fgsA, ' FGS B =', fgsB READ (1, '(A)'), firstline1 READ (1, '(A)'), firstline1 END IF dataline1 = 5 dataline2 = 5 ! Prepare to maintain line positions first = .TRUE. seconds2 = -999 ! Initialise to deliberately unreasonable value TYPE *, ' Processing. This may take several minutes' 1 IF (fine) THEN READ (1, '(A20, 2X, 6F16.5)', IOSTAT=stat1, END=4), X time, rpmtxa1, rpmtxb1, rpmtya1, rpmtyb1, theta_a1, theta_b1 dataline1 = dataline1+1 IF (stat1 .EQ. 0) THEN ! First interpret time string READ (time, '(1X,I2,1X,I3,1X,I2,1X,I2,1X,F6.3)', IOSTAT=stat3) year, day, hour, min ,sec IF (stat3 .NE. 0) THEN TYPE *, ' ERROR IN DATE STRING AT LINE ', dataline1 GO TO 1 END IF IF (first) day0 = day first = .FALSE. hours = (sec/60.0+min)/60.0+hour + 24.0*(day-day0) seconds = 3600*hours IF (seconds-seconds2 .GT. 0.51) GO TO 2 IF (seconds2-seconds .GT. 0.51) GO TO 1 GO TO 5 ! Times are synchronised - Produce output ELSE ! Bad stat1 TYPE *, ' BAD DATA AT LINE ', dataline1, ' IN ', fes1, '...IGNORED' GO TO 1 END IF 2 READ (2, '(A20, 2X, 6F16.5)', IOSTAT= stat2, END=4), X time2, rpmtxa2, rpmtxb2, rpmtya2, rpmtyb2, theta_a2, theta_b2 dataline2 = dataline2+1 IF (stat2 .EQ. 0) THEN READ (time2, '(1X,I2,1X,I3,1X,I2,1X,I2,1X,F6.3)', IOSTAT=stat3) year, day, hour, min ,sec IF (stat3 .NE. 0) THEN TYPE *, ' ERROR IN DATE STRING AT LINE ', dataline2 GO TO 2 END IF IF (first) day0 = day first = .FALSE. hours2 = (sec/60.0+min)/60.0+hour + 24.0*(day-day0) seconds2 = 3600*hours2 IF (seconds-seconds2 .GT. 0.51) GO TO 2 IF (seconds2-seconds .GT. 0.51) GO TO 1 ELSE ! Bad stat2 TYPE *, ' BAD DATA AT LINE ', dataline2, ' IN ', fes2,'...IGNORED' GO TO 2 END IF C Produce a line of output. 5 pmtxa1 = INT(rpmtxa1) pmtxb1 = INT(rpmtxb1) pmtya1 = INT(rpmtya1) pmtyb1 = INT(rpmtyb1) Call FGS_TO_VEH (fgsA, theta_a1, theta_b1, pmtxa1, pmtya1, pmtxb1, pmtyb1, X gm, alpha, beta, tvs, X alever(fgsA), blever(fgsA), aoffset(fgsA), boffset(fgsA), k1x(fgsA), k1y(fgsA), X fgs_ima_vec, fgs_obj_vec, fgs_pos, fgs_veh_vec1) pmtxa2 = INT(rpmtxa2) pmtxb2 = INT(rpmtxb2) pmtya2 = INT(rpmtya2) pmtyb2 = INT(rpmtyb2) Call FGS_TO_VEH (fgsB, theta_a2, theta_b2, pmtxa2, pmtya2, pmtxb2, pmtyb2, X gm, alpha, beta, tvs, X alever(fgsB), blever(fgsB), aoffset(fgsB), boffset(fgsB),k1x(fgsB), k1y(fgsB), X fgs_ima_vec, fgs_obj_vec, fgs_pos, fgs_veh_vec2) WRITE (3, '(1X, F8.5, 4F12.3, 2X, A22)') hours, fgs_veh_vec1(2), X fgs_veh_vec1(3), fgs_veh_vec2(2), fgs_veh_vec2(3), time IF (MOD(dataline1,500) .EQ. 0) TYPE *, ' Done ', dataline1, ' lines.' ELSE ! Coarse 3 READ (1, '(A20, 2X, 4F16.5)', IOSTAT=stat1, END=4) X time, theta_a1, theta_b1, theta_a2, theta_b2 dataline1 = dataline1 + 1 IF (MOD(dataline1,500) .EQ. 0) TYPE *, ' Done ', dataline1, ' lines.' IF (stat1 .EQ. 0) THEN Call FGS_TO_VEH (fgsA, theta_a1,theta_b1, 0.0, 0.0, 0.0, 0.0, gm, alpha, beta, tvs, X alever(fgsA), blever(fgsA), aoffset(fgsA), boffset(fgsA), k1x(fgsA), k1y(fgsA), X fgs_ima_vec, fgs_obj_vec, fgs_pos, fgs_veh_vec1) Call FGS_TO_VEH (fgsB, theta_a2,theta_b2, 0.0, 0.0, 0.0, 0.0, gm, alpha, beta, tvs, X alever(fgsB), blever(fgsB), aoffset(fgsB), boffset(fgsB), k1x(fgsB), k1y(fgsB), X fgs_ima_vec, fgs_obj_vec, fgs_pos, fgs_veh_vec2) READ (time, '(1X,I2,1X,I3,1X,I2,1X,I2,1X,F6.3)', IOSTAT=stat3) year, day, hour, min ,sec IF (stat3 .EQ. 0) THEN ! Allow for 24-hour wraparound IF (first) day0 = day first = .FALSE. hours = (sec/60.0+min)/60.0+hour + 24.0*(day-day0) WRITE (3, '(1X, F8.5, 4F12.3, 2X, A22)') hours, fgs_veh_vec1(2), X fgs_veh_vec1(3), fgs_veh_vec2(2), fgs_veh_vec2(3), time dataline1 = dataline1+1 IF (MOD(dataline1,500) .EQ. 0) TYPE *, ' Done ', dataline1, ' lines.' ELSE TYPE *, ' ERROR IN DATE STRING AT LINE ', dataline1 GO TO 3 END IF ELSE ! Bad stat1 TYPE *, ' BAD DATA AT LINE ', dataline1, ' IN ', guide, '...IGNORED' GO TO 1 END IF END IF ! fine or coarse GO TO 1 4 CLOSE (2) CLOSE(3) RETURN END C***************************************************************************** SUBROUTINE READSCHF(schf, alever, blever,tvs, k1x, k1y, aoffset, boffset, gm, X alpha, beta, found) C***************************************************************************** C C Read the mini-SCHF file and set up variables C IMPLICIT NONE INTEGER i, j, f, schf CHARACTER*(20) quote CHARACTER*80 text, datatext REAL*8 alever(3), blever(3), tvs(3,3,3), k1x(3), k1y(3) REAL*8 k3x, k3y, aoffset(3), boffset(3) REAL*8 gm(3), alpha(6,6,3), beta(6,6,3) LOGICAL found quote = 'K-FACTORS' Call LOOKFOR (schf, text, quote, found) IF (.NOT. found) GO TO 999 READ (schf, *) k1x(1), k3x, k1y(1), k3y READ (schf, *) k1x(2), k3x, k1y(2), k3y READ (schf, *) k1x(3), k3x, k1y(3), k3y quote = 'OPTICAL MAGNIFICATION' Call LOOKFOR (schf, text, quote, found) IF (.NOT. found) GO TO 999 READ (schf, *) gm quote = 'OPTCOF' Call LOOKFOR (schf, text, quote, found) IF (.NOT. found) GO TO 999 C Undefined alpha and beta values are set to zero DO i = 1, 6 DO j = 1, 6 DO f = 1, 3 alpha(i,j,f) = 0.0D0 beta(i,j,f) = 0.0D0 END DO END DO END DO DO f = 1, 3 READ (schf, *) alpha(2,1,f), alpha(2,2,f), alpha(3,1,f) READ (schf, *) alpha(1,3,f), alpha(4,1,f), alpha(2,3,f) READ (schf, *) alpha(6,1,f), alpha(4,3,f), alpha(2,5,f) READ (schf, *) alpha(1,4,f), alpha(3,2,f) READ (schf, *) beta(1,2,f), beta(2,2,f), beta(1,3,f) READ (schf, *) beta(3,1,f), beta(1,4,f), beta(3,2,f) READ (schf, *) beta(1,6,f), beta(3,4,f), beta(5,2,f) READ (schf, *) beta(4,1,f), beta(2,3,f) END DO quote = 'TVS MATRICES' Call LOOKFOR (schf, text, quote, found) IF (.NOT. found) GO TO 999 READ (schf, *) tvs quote = 'ALEVER' Call LOOKFOR (schf, text, 'ALEVER', found) IF (.NOT. found) GO TO 999 READ (schf, *) alever quote = 'BLEVER' Call LOOKFOR (schf, text, quote, found) IF (.NOT. found) GO TO 999 READ (schf, *) blever quote = 'STAR SELECTOR OFFSETS' Call LOOKFOR(schf, text, quote, found) IF (.NOT. found) GO TO 999 READ (schf, '(A)') text ! Contains SSOFFA= datatext = text(8:80) ! Trim this off READ(datatext,*) aoffset READ (schf, '(A)') text ! Contains SSOFFB= datatext = text(8:80) ! Trim this off READ(datatext,*) boffset RETURN 999 TYPE *, quote, ' NOT FOUND ' ! found is false RETURN END C*************************************************************************** SUBROUTINE LOOKFOR( file, textline, quote, found) C*************************************************************************** IMPLICIT NONE INTEGER*4 file CHARACTER*80 textline CHARACTER*(*) quote LOGICAL found C C Finds and returns a line of text containing the quoted string C found = .false. DO WHILE (.NOT. found) READ(file, '(A)', END=1) textline found = (INDEX(textline, quote) .GT. 0) END DO 1 RETURN END C**************************************************************************** SUBROUTINE TARGETPOS C**************************************************************************** C C Locate Target in V2V3 C Using Guide star RA and DECs plus V2V3s find telescope attitude C Then convert Target RA and DEC to V2V3 or vice-versa IMPLICIT NONE REAL*8 raP, decP, raS, decS, v2P, v2S, v3P, v3S, v2T, v3T REAL*8 v1ra, v1dec, roll, raT, decT, xc, yc, zc, vx, vy, vz, v2, v3 REAL*8 velx, vely, velz,hour,oldhour,oldvelx,oldvely,oldvelz REAL*8 newhour, newvelx, newvely, newvelz REAL*8 fhour, sec, craP, cdecP, craS, cdecS, rms INTEGER*4 yr, day, day0, hr, min, p, velerr, tline CHARACTER*40 starfile, velocities, results, timestring, option, wantvel,targ CHARACTER*80 targstring CHARACTER*80 TARGNAME LOGICAL found, filepresent, velflag CHARACTER*1 OPT C Read in guide star information 29 TYPE *, ' Give file name containing star positions' READ '(A)', starfile INQUIRE(file=starfile,exist=filepresent) IF (filepresent) then OPEN (1, file= starfile,type='old') READ (1,'(A)') option READ (1,*) raP, decP READ (1,*) raS, decS p=0 c Open output file 55 p=p+1 READ (1,'(A)',end=20) targstring IF ((INDEX(option,'v2').gt.0).or.(INDEX(option,'V2').gt.0)) THEN opt='v' READ (targstring,*) raT, decT targ=TARGNAME(targstring) results= targ(1:INDEX(targ, ' ') -1) // '.v2v3' OPEN (4, file = results) write (4,*) ' Telescope Attitude (V1 axis) Target Guidestar sep (o-c)' write (4,*) ' Time (hours) RA DEC Roll V2 V3 (diff in arcsec) year/day/h/m/s' write (4,*) ELSE IF ((INDEX(option,'ra').gt.0).or.(INDEX(option,'RA').gt.0)) THEN opt='r' READ (targstring,*) v2T, v3T targ=TARGNAME(targstring) results= targ(1:INDEX(targ, ' ') -1) // '.radec' OPEN (4, file = results) write (4,*) ' Telescope Attitude (V1 axis) Target Guidestar sep (o-c)' write (4,*) ' Time (hours) RA DEC Roll RA DEC (diff in arcsec)' write (4,*) ELSE TYPE *,' THE FIRST LINE OF THE STAR FILE MUST CONTAIN EITHER' TYPE *,' AN UPPER OR LOWER CASE RA OR V2 TO INDICATE WHICH WAY TO' TYPE *,' TRANSFORM RESULT. PLEASE EDIT THAT FILE AND RERUN.' STOP ENDIF ELSE type *,' THERE IS NO FILE NAMED ',starfile type *,' PLEASE RE-ENTER' goto 29 ENDIF C Open velocities files IF (p.eq.1) then TYPE *, 'Do you wish to perform velocity aberration correction?' READ '(A)', wantvel velflag = ((wantvel(1:1) .EQ. 'y') .OR. (wantvel(1:1) .EQ. 'Y') ) IF (velflag) THEN 31 TYPE * , 'Give the velocity file' READ '(A)', velocities INQUIRE(file=velocities,exist=filepresent) IF (.not.filepresent) then type *,' THERE IS NO FILE NAMED ',velocities type *,' PLEASE RE-ENTER' goto 31 ENDIF ELSE ENDIF ENDIF IF (velflag) THEN OPEN (2,file=velocities,type='old') ENDIF OPEN (3,file='guidestars.v23',type='old') IF (velflag) THEN Read (2,900,end=19, iostat=velerr) yr,day,hr,min,sec,oldvelx,oldvely,oldvelz day0 = day ! Allow for 24-hour wraparound oldhour = (sec/60+min)/60+hr + 24.0*(day-day0) READ (2,100, end=9, err=19) yr, day, hr, min, sec, * newvelx, newvely, newvelz 900 FORMAT (///// 1X,I2, 1X, I3, 2(1X, I2), 1X, F6.3, 3x, 3(f16.5)) 100 FORMAT (1X,I2, 1X, I3, 2(1X, I2), 1X, F6.3, 3x, 3(f16.5)) newhour = (sec/60+min)/60+hr + 24.0*(day-day0) END IF tline = 0 10 READ (3,23,end=9) hour, v2p, v3p, v2s, v3s, timestring 23 format (f9.5, 4(f12.3), A22) IF (velflag) THEN found = .false. DO WHILE (.NOT.found) IF (newhour .GE. hour) THEN found = .true. fhour = (hour-oldhour)/(newhour-oldhour) velx = oldvelx + fhour*(newvelx-oldvelx) vely = oldvely + fhour*(newvely-oldvely) velz = oldvelz + fhour*(newvelz-oldvelz) ELSE oldvelx = newvelx oldvely = newvely oldvelz = newvelz oldhour = newhour 32 READ (2,100, end=21, iostat=velerr) yr, day, hr, min, sec, * newvelx, newvely, newvelz IF (velerr .NE. 0) THEN TYPE *, ' ERROR READING VELOCITY FILE - INTERPOLATING TO ', timestring GO TO 32 END IF newhour = (sec/60+min)/60+hr + 24.0*(day-day0) END IF END DO END IF C Apply velocity aberration correction Call RADEC_U (raP, decP, xc, yc, zc) IF (velflag) Call NONREL(xc, yc, zc, velx, vely, velz) Call U_RADEC(xc, yc, zc, craP, cdecP) Call RADEC_U (raS, decS, xc, yc, zc) IF (velflag) Call NONREL(xc, yc, zc, velx, vely, velz) Call U_RADEC(xc, yc, zc, craS, cdecS) Call ATTITUDE (timestring,craP, cdecP, v2P, v3P, craS, cdecS, v2S, v3S, * v1ra, v1dec, roll, rms) C FIND TARGET LOCATION IF (opt.eq.'v') THEN C Compute V2,V3 of Target Call RADEC_U (raT,decT, xc, yc, zc) IF (velflag) Call NONREL(xc, yc, zc, velx, vely, velz) Call CEL_TEL (v1ra, v1dec, roll, xc, yc, zc, vx, vy, vz) Call U_V2V3 (vx, vy, vz, v2, v3) IF (v2**2+v3**2 .GT. 360000) TYPE *, ' THE TARGET IS OUTSIDE HST FOV' C Put output in file "target.v2v3" WRITE (4,'(1X, F12.8, 3F12.6, 2F10.3, F12.3, A22)') * hour, v1ra, v1dec, roll, v2, v3, rms, timestring ELSE IF (opt.eq.'r') THEN C Compute heliocentric RA, and DEC of Target Call V2V3_U (v2T, v3T, vx, vy, vz) Call TEL_CEL (v1ra, v1dec, roll, vx, vy, vz, xc, yc, zc) IF (velflag) Call NONREL(xc, yc, zc, -velx, -vely, -velz) Call U_RADEC (xc, yc, zc, raT, decT) C Put output in file "target.radec" WRITE (4,'(1X, F12.8, 3F15.8, 2X, 2F15.8, F12.3, A22)') * hour, v1ra, v1dec, roll, raT, decT, rms, timestring ENDIF tline = tline+1 IF (MOD(tline,500) .EQ. 0) TYPE *, ' Done ', tline, ' target positions.' GOTO 10 21 TYPE *, ' REACHED END OF VELOCITY FILE' 9 IF (velflag) THEN CLOSE(2) ENDIF CLOSE(3) Type *,'Results written to the file ',results 19 GOTO 55 20 CLOSE(1) RETURN END C************************************************************************** SUBROUTINE RADEC_U (ra, dec, vx, vy, vz) C************************************************************************** C C Express RA and DEC direction as a unit vector. RA and DEC in degrees IMPLICIT NONE REAL*8 ra, dec, vx, vy, vz vx = DCOSD(ra)*DCOSD(dec) vy = DSIND(ra)*DCOSD(dec) vz = DSIND(dec) RETURN END C*************************************************************************** SUBROUTINE U_RADEC (xc, yc, zc, ra, dec) C*************************************************************************** IMPLICIT NONE REAL*8 xc, yc, zc, ra, dec C Convert celestial coordinate unit vector to RA and DEC in degrees ra = DATAN2D(yc, xc) IF (ra .LT. 0.0) ra = ra+360.0 dec = DASIND(zc) RETURN END C************************************************************************** SUBROUTINE V2V3_U (v2, v3, vx, vy, vz) C************************************************************************** C C Express V2, V3 as unit vector in telescope space. V2V3 in arcsec IMPLICIT NONE REAL*8 v2, v3, vx, vy, vz REAL*8 rho rho = DCOSD(v3/3600.0) vx = rho*DCOSD(v2/3600) vy = rho*DSIND(v2/3600) vz = DSIND(v3/3600) RETURN END C************************************************************************** SUBROUTINE U_V2V3 (vx, vy, vz, v2, v3) C************************************************************************** C C Convert unit vector in telescope space to V2 and V3 in arcsec IMPLICIT NONE REAL*8 vx, vy, vz, v2, v3 v2 = 3600*DATAN2D(vy, vx) v3 = 3600*DASIND(vz) RETURN END C************************************************************************* SUBROUTINE CEL_TEL (v1ra, v1dec, roll, xc, yc, zc, vx, vy, vz) C************************************************************************* C C Transform from unit vector in celestial coords to telescope coords C This is a unit vector in the V1,V2,V3 system REAL*8 v1ra, v1dec, roll, xc, yc, zc, vx, vy, vz REAL*8 x1, y1, z2 x1 = xc*DCOSD(v1ra) + yc*DSIND(v1ra) y1 = -xc*DSIND(v1ra) + yc*DCOSD(v1ra) vx = x1*DCOSD(v1dec) + zc*DSIND(v1dec) z2 = -x1*DSIND(v1dec) + zc*DCOSD(v1dec) vy = y1*DCOSD(roll) - z2*DSIND(roll) vz = y1*DSIND(roll) + z2*DCOSD(roll) RETURN END C************************************************************************* SUBROUTINE TEL_CEL (v1ra, v1dec, roll, vx, vy, vz, xc, yc, zc) C************************************************************************* C C Transform from unit vector in v1,v2,v3 to celestial xc, yc, zc REAL*8 v1ra, v1dec, roll, vx, vy, vz, xc, yc, zc REAL*8 y2, z2, x1 y2 = vy*DCOSD(roll) + vz*DSIND(roll) z2 = -vy*DSIND(roll) + vz*DCOSD(roll) x1 = vx*DCOSD(v1dec) - z2*DSIND(v1dec) zc = vx*DSIND(v1dec) + z2*DCOSD(v1dec) xc = x1*DCOSD(v1ra) - y2*DSIND(v1ra) yc = x1*DSIND(v1ra) + y2*DCOSD(v1ra) RETURN END C************************************************************************* SUBROUTINE NONREL (ax, ay, az, vx, vy, vz) C************************************************************************* C C Corrects for velocity aberration of spacecraft motion C C (ax, ay, az) is true direction as a unit vector C (vx, vy, vz) is velocity vector in metres/sec C Apparent viewing direction returned in (ax, ay, az) C IMPLICIT NONE REAL*8 ax, ay, az, vx, vy, vz, a, c PARAMETER (c=2.9979D8) ! Velocity of light in m/s C Apply velocity aberration correction ax = ax + vx/c ay = ay + vy/c az = az + vz/c a = DSQRT(ax**2+ay**2+az**2) ax = ax/a ay = ay/a az = az/a RETURN END C************************************************************************** SUBROUTINE ATTITUDE (timestring, ra, da, v2a, v3a, rb, db, v2b, v3b, * v1ra, v1dec, roll, rms) C************************************************************************** C The guide star calculation - Given RA and DEC of two stars and C their V2 V3 positions, find telescope attitude, v1ra, v1dec, roll C V2 and V3 are always given in arcsec. Other angles in degrees. C C New method using least-square deviation IMPLICIT NONE C Common block description REAL*8 ra, da, v2a, v3a, x0a, y0a, z0a, x1a, y1a, z1a REAL*8 x2a, y2a, z2a, x3a, y3a, z3a, xva, yva, zva REAL*8 rb, db, v2b, v3b, x0b, y0b, z0b, x1b, y1b, z1b REAL*8 x2b, y2b, z2b, x3b, y3b, z3b, xvb, yvb, zvb COMMON /STARPOS/x0a, y0a, z0a, x1a, y1a, z1a, x2a, y2a, z2a, 2 x3a, y3a, z3a, xva, yva, zva, 3 x0b, y0b, z0b, x1b, y1b, z1b, x2b, y2b, z2b, 4 x3b, y3b, z3b, xvb, yvb, zvb CHARACTER*40 timestring REAL*8 v1ra, v1dec, roll REAL*8 angles(3), ftol, fret, rad2arcs, rms INTEGER*4 nvec, iter PARAMETER (rad2arcs = 2.062D5) ! Radians to arcsec conversion C First get approximate solution Call APPROX (timestring, ra, da, rb, db, v2a, v3a, v2b, v3b, X v1ra, v1dec, roll) C Now refine the solution by minimizing the sum of separations squared C between the calculated and actual direction vectors C using Numerical Recipes routine FRPRMN (Fletcher-Reeves-Polak-Ribiere) C Unit vectors defined by V2V3 values Call V2V3_U(v2a, v3a, xva, yva, zva) Call V2V3_U(v2b, v3b, xvb, yvb, zvb) angles(1) = v1ra angles(2) = v1dec angles(3) = roll C Unit vector of stars in GCI Call RADEC_U (ra, da, x0a, y0a, z0a) Call RADEC_U (rb, db, x0b, y0b, z0b) nvec = 3 ftol = 1.0D-8 ! Should try for 1/10 mas Call FRPRMN (angles, nvec, ftol, iter, fret) rms = DSQRT(fret/2)*rad2arcs IF (rms .GT. 2.0) TYPE 100, timestring, rms 100 FORMAT (A40,'BAD FIT. RMS DEVIATION ',F8.3,' arcsec'/) v1ra = angles(1) v1dec = angles(2) roll = angles(3) rms = DSQRT(fret/2)*rad2arcs RETURN END C---------------------------------------------------------------------- REAL*8 FUNCTION FUNC(angles) C Separation-squared to be minimized IMPLICIT NONE C Common block description REAL*8 x0a, y0a, z0a, x1a, y1a, z1a REAL*8 x2a, y2a, z2a, x3a, y3a, z3a, xva, yva, zva REAL*8 x0b, y0b, z0b, x1b, y1b, z1b REAL*8 x2b, y2b, z2b, x3b, y3b, z3b, xvb, yvb, zvb COMMON /STARPOS/x0a, y0a, z0a, x1a, y1a, z1a, x2a, y2a, z2a, 2 x3a, y3a, z3a, xva, yva, zva, 3 x0b, y0b, z0b, x1b, y1b, z1b, x2b, y2b, z2b, 4 x3b, y3b, z3b, xvb, yvb, zvb REAL*8 angles(3) REAL*8 cph, sph, cth, sth, cps, sps, dsa, dsb C Collect frequently used sines and cosines cph = COSD(angles(1)) sph = SIND(angles(1)) cth = COSD(angles(2)) sth = SIND(angles(2)) cps = COSD(angles(3)) sps = SIND(angles(3)) C Perform transformation to telescope coordinates x1a = x0a*cph + y0a*sph y1a = -x0a*sph + y0a*cph z1a = z0a x2a = x1a*cth + z1a*sth y2a = y1a z2a = -x1a*sth + z1a*cth x3a = x2a y3a = y2a*cps - z2a*sps z3a = y2a*sps + z2a*cps C Same for star-b x1b = x0b*cph + y0b*sph y1b = -x0b*sph + y0b*cph z1b = z0b x2b = x1b*cth + z1b*sth y2b = y1b z2b = -x1b*sth + z1b*cth x3b = x2b y3b = y2b*cps - z2b*sps z3b = y2b*sps + z2b*cps C Square of length between vectors dsa = (x3a-xva)**2 + (y3a-yva)**2 + (z3a-zva)**2 dsb = (x3b-xvb)**2 + (y3b-yvb)**2 + (z3b-zvb)**2 func = dsa + dsb RETURN END C------------------------------------------------------------------------- SUBROUTINE DFUNC (angles, df) C Provides differentials of the separation-squared function IMPLICIT NONE REAL*8 angles(3), df(3), pi REAL*8 cph, sph, cth, sth, cps, sps REAL*8 dx3adph, dy3adph, dz3adph REAL*8 dx3adth, dy3adth, dz3adth REAL*8 dx3adps, dy3adps, dz3adps REAL*8 dx3bdph, dy3bdph, dz3bdph REAL*8 dx3bdth, dy3bdth, dz3bdth REAL*8 dx3bdps, dy3bdps, dz3bdps REAL*8 dsxa, dsya, dsza, dsxb, dsyb, dszb REAL*8 dsadph, dsadth, dsadps, dsbdph, dsbdth, dsbdps REAL*8 x0a, y0a, z0a, x1a, y1a, z1a REAL*8 x2a, y2a, z2a, x3a, y3a, z3a, xva, yva, zva REAL*8 x0b, y0b, z0b, x1b, y1b, z1b REAL*8 x2b, y2b, z2b, x3b, y3b, z3b, xvb, yvb, zvb COMMON /STARPOS/x0a, y0a, z0a, x1a, y1a, z1a, x2a, y2a, z2a, 2 x3a, y3a, z3a, xva, yva, zva, 3 x0b, y0b, z0b, x1b, y1b, z1b, x2b, y2b, z2b, 4 x3b, y3b, z3b, xvb, yvb, zvb pi = 4.0*DATAN(1.0D0) C Collect frequently used sines and cosines cph = COSD(angles(1)) sph = SIND(angles(1)) cth = COSD(angles(2)) sth = SIND(angles(2)) cps = COSD(angles(3)) sps = SIND(angles(3)) C Calculate partial differentials dx3adph = y1a*cth dy3adph = -x1a*cps + y1a*sps*sth dz3adph = -x1a*sps - y1a*cps*sth dx3adth = z2a dy3adth = x2a*sps dz3adth = -x2a*cps dx3adps = 0.0 dy3adps = -z3a dz3adps = y3a C b contribution dx3bdph = y1b*cth dy3bdph = -x1b*cps + y1b*sps*sth dz3bdph = -x1b*sps - y1b*cps*sth dx3bdth = z2b dy3bdth = x2b*sps dz3bdth = -x2b*cps dx3bdps = 0.0 dy3bdps = -z3b dz3bdps = y3b C Vector differences dsxa = x3a - xva dsya = y3a - yva dsza = z3a - zva dsxb = x3b - xvb dsyb = y3b - yvb dszb = z3b - zvb C Now build the total differentials dsadph = dsxa*dx3adph + dsya*dy3adph + dsza*dz3adph dsadth = dsxa*dx3adth + dsya*dy3adth + dsza*dz3adth dsadps = dsxa*dx3adps + dsya*dy3adps + dsza*dz3adps dsbdph = dsxb*dx3bdph + dsyb*dy3bdph + dszb*dz3bdph dsbdth = dsxb*dx3bdth + dsyb*dy3bdth + dszb*dz3bdth dsbdps = dsxb*dx3bdps + dsyb*dy3bdps + dszb*dz3bdps df(1) = pi/90*(dsadph + dsbdph) df(2) = pi/90*(dsadth + dsbdth) df(3) = pi/90*(dsadps + dsbdps) RETURN END C*********************************************************************** SUBROUTINE APPROX (timestring, ra, da, rb, db, v2a, v3a, v2b, v3b, X v1ra, v1dec, roll) C*********************************************************************** C Get approximate attitude to start accurate fit IMPLICIT NONE CHARACTER*40 timestring REAL*8 ra, da, rb, db, v2a, v3a, v2b, v3b, v1ra, v1dec, roll REAL*8 dec, de, dn, dv2, dv3, ea, na REAL*8 sep1, sep2, scale dec = 0.5*(da+db) ! Mean declination de = 3600*(rb-ra)*COSD(dec) ! In arcsec dn = 3600*(db-da) dv2 = v2b-v2a dv3 = v3b-v3a sep1 = SQRT(de**2+dn**2) sep2 = SQRT(dv2**2 +dv3**2) scale = sep1/sep2 ! Allow for small scale difference IF ((scale .LT. 0.99) .OR. (scale .GT. 1.01)) * TYPE 100, timestring, scale 100 FORMAT (A40, 'ABNORMAL SCALE ',F8.3) roll = ATAN2D(de*dv3-dn*dv2, de*dv2+dn*dv3 ) ea = (v2a*COSD(roll) + v3a*SIND(roll))/3600 ! Back to degree\s na = (-v2a*SIND(roll) + v3a*COSD(roll))/3600 v1ra = ra-scale*ea/COSD(dec) v1dec = da - scale*na RETURN END C*********************************************************************** SUBROUTINE frprmn(p,n,ftol,iter,fret) C*********************************************************************** INTEGER iter,n,NMAX,ITMAX REAL*8 func,fret,ftol,p(n),EPS PARAMETER (NMAX=50,ITMAX=200,EPS=1.e-10) CU USES dfunc,func,linmin INTEGER its,j REAL*8 dgg,fp,gam,gg,g(NMAX),h(NMAX),xi(NMAX) fp=func(p) call dfunc(p,xi) do j=1,n g(j)=-xi(j) h(j)=g(j) xi(j)=h(j) enddo do its=1,ITMAX iter=its call linmin(p,xi,n,fret) if(2.*abs(fret-fp).le.ftol*(abs(fret)+abs(fp)+EPS))return fp=func(p) call dfunc(p,xi) gg=0. dgg=0. do j=1,n gg=gg+g(j)**2 C dgg=dgg+xi(j)**2 dgg=dgg+(xi(j)+g(j))*xi(j) enddo if(gg.eq.0.)return gam=dgg/gg do j=1,n g(j)=-xi(j) h(j)=g(j)+gam*h(j) xi(j)=h(j) enddo enddo pause 'frpr maximum iterations exceeded' return END C*************************************************************************** SUBROUTINE linmin(p,xi,n,fret) C*************************************************************************** IMPLICIT REAL*8 (A-H, O-Z) INTEGER n,NMAX REAL*8 fret,p(n),xi(n),TOL COMMON /f1com/ ncom,pcom,xicom PARAMETER (NMAX=50,TOL=1.e-4) CU USES brent,f1dim,mnbrak INTEGER j,ncom REAL*8 pcom(NMAX),xicom(NMAX),ax,bx,fa,fb,fx,xmin,xx,brent external f1dim ncom=n do j=1,n pcom(j)=p(j) xicom(j)=xi(j) enddo ax=0. xx=1. call mnbrak(ax,xx,bx,fa,fx,fb,f1dim) fret=brent(ax,xx,bx,f1dim,TOL,xmin) do j=1,n xi(j)=xmin*xi(j) p(j)=p(j)+xi(j) enddo return END C----------------------------------------------------------------------------- FUNCTION brent(ax,bx,cx,f,tol,xmin) IMPLICIT REAL*8 (A-H, O-Z) INTEGER ITMAX REAL*8 brent,ax,bx,cx,tol,xmin,f,CGOLD,ZEPS EXTERNAL f PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0e-10) INTEGER iter REAL*8 a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm a=min(ax,cx) b=max(ax,cx) v=bx w=v x=v e=0. fx=f(x) fv=fx fw=fx do iter=1,ITMAX xm=0.5*(a+b) tol1=tol*abs(x)+ZEPS tol2=2.*tol1 if(abs(x-xm).le.(tol2-.5*(b-a))) goto 3 if(abs(e).gt.tol1) then r=(x-w)*(fx-fv) q=(x-v)*(fx-fw) p=(x-v)*q-(x-w)*r q=2.*(q-r) if(q.gt.0.) p=-p q=abs(q) etemp=e e=d if(abs(p).ge.abs(.5*q*etemp).or.p.le.q*(a-x).or. * p.ge.q*(b-x)) goto 1 d=p/q u=x+d if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x) goto 2 endif 1 if(x.ge.xm) then e=a-x else e=b-x endif d=CGOLD*e 2 if(abs(d).ge.tol1) then u=x+d else u=x+sign(tol1,d) endif fu=f(u) if(fu.le.fx) then if(u.ge.x) then a=x else b=x endif v=w fv=fw w=x fw=fx x=u fx=fu else if(u.lt.x) then a=u else b=u endif if(fu.le.fw .or. w.eq.x) then v=w fv=fw w=u fw=fu else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then v=u fv=fu endif endif enddo pause 'brent exceed maximum iterations.' 3 xmin=x brent=fx return END C-------------------------------------------------------------------------- FUNCTION f1dim(x) IMPLICIT REAL*8 (A-H, O-Z) INTEGER NMAX REAL*8 f1dim,func,x COMMON /f1com/ ncom,pcom,xicom PARAMETER (NMAX=50) CU USES func INTEGER j,ncom REAL*8 pcom(NMAX),xicom(NMAX),xt(NMAX) do j=1,ncom xt(j)=pcom(j)+x*xicom(j) enddo f1dim=func(xt) return END c********************************************************************** CHARACTER*(*) FUNCTION TARGNAME (targstring) IMPLICIT NONE INTEGER l, first, last CHARACTER*1 c, tab CHARACTER*(*) targstring DATA tab /9/ l = LEN(targstring) c = targstring(l:l) DO WHILE (c .EQ. ' ') l = l-1 c = targstring(l:l) END DO last = l ! End of target name DO WHILE (( c .NE. ' ') .AND. (c .NE. tab)) l = l-1 c = targstring(l:l) END DO first = l+1 ! Beginning of target name targname = targstring(first:last) RETURN END C************************************************************************** SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func) C************************************************************************** IMPLICIT REAL*8 (A-H, O-Z) REAL*8 ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY EXTERNAL func PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.e-20) REAL*8 dum,fu,q,r,u,ulim fa=func(ax) fb=func(bx) if(fb.gt.fa)then dum=ax ax=bx bx=dum dum=fb fb=fa fa=dum endif cx=bx+GOLD*(bx-ax) fc=func(cx) 1 if(fb.ge.fc)then r=(bx-ax)*(fb-fc) q=(bx-cx)*(fb-fa) u=bx-((bx-cx)*q-(bx-ax)*r)/(2.*sign(max(abs(q-r),TINY),q-r)) ulim=bx+GLIMIT*(cx-bx) if((bx-u)*(u-cx).gt.0.)then fu=func(u) if(fu.lt.fc)then ax=bx fa=fb bx=u fb=fu return else if(fu.gt.fb)then cx=u fc=fu return endif u=cx+GOLD*(cx-bx) fu=func(u) else if((cx-u)*(u-ulim).gt.0.)then fu=func(u) if(fu.lt.fc)then bx=cx cx=u u=cx+GOLD*(cx-bx) fb=fc fc=fu fu=func(u) endif else if((u-ulim)*(ulim-cx).ge.0.)then u=ulim fu=func(u) else u=cx+GOLD*(cx-bx) fu=func(u) endif ax=bx bx=cx cx=u fa=fb fb=fc fc=fu go to 1 endif return END c**************************************************************************** SUBROUTINE FGS_TO_VEH (FGS_BAY,THETA_A, THETA_B, PMTXA, PMTYA, & PMTXB, PMTYB,GM,ALPHA,BETA,TVS, & ALEVER, BLEVER, aoffst, boffst, K1X, K1Y, & FGS_IMA_VEC, FGS_OBJ_VEC,FGS_POS,FGS_VEH_VEC) c**************************************************************************** C -Form guide star position unit vector in vehicle space. C C C------------------------------------------------------------------------------C C C c c Purpose: Form the instantaneous guide star position unit vector c in Vehicle Space (when FGS is in fine-lock). C c Input data: PMT counts, FGS encoder angles C FGS constants: (to be obtained from PDB) C FGS lever arms and encoder offsets C FGS K factors C FGS distortion coefficients C FGS to Vehicle frame transformation matrix C c Output data: FGS fine-lock vector in V1V2V3 space (unit vector) C C C ARGUMENT LIST: C C ARGUMENT NAME TYPE USE DESCRIPTION C --------------- ---- --- ------------------------------ c FGS_bay I FGS bay c FGS constants: c ------------- C ALEVER R*8 I FGS star selector a deflection c angle (arcminutes) C BLEVER R*8 I FGS star selector b deflection c angle (arcminutes) C AOFFST R*8 I FGS star selector a offset angle C (arcminutes) C BOFFST R*8 I FGS star selector b offset angle c (arcminutes) C K1X R*8 I FGS interferometer c coefficient (arc-min) C K1Y R*8 I FGS interferometer c coefficient (arc-min) C K0X R*8 I FGS interferometer c coefficient (arc-min) C K0Y R*8 I FGS interferometer c coefficient (arc-min) C ALPHA(NX,NY,FGS_BAY) R*8 I FGS x-axis optical distortion c coefficients: c index 1 - x coefficient c index 2 - y coefficient c C BETA(MX,MY,FGS_BAY) R*8 I FGS y-axis optical distortion c coefficients: c index 1 - x coefficient c index 2 - y coefficient c C NX I*4 I Mumber of x coefficient terms for c x component c C NY I*4 I Number of y coefficient terms for c x component c C MX I*4 I Number of x coefficient terms for c y component c C MY I*4 I Number of y coefficient terms for c y component C GM(3) R*8 I FGS magnification factor (FGS_BAY) c C TVS(3,3,3) R*8 I FGS to ST vehicle transformation c array c (row,column,bay) c TVS_MAT(3,3) R*8 I the corresponding matrix for c the FGS under study C C FGS measurement data: c -------------------- C C PMTXA I*4 I Photomultiplier tube (PMT) c A x counts C PMTYA I*4 I PMT A Y counts C PMTXB I*4 I PMT B X counts C PMTYB I*4 I PMT B Y counts C C THETA_A R*8 I Star selector angle A (arcseconds) C THETA_B R*8 I Star selector angle B (arcseconds) c Intermediate and final calculation output: c ------------------------------------------ C C RX R*8 O X Fine error signal C RY R*8 O Y Fine error signal C C FGS_POS(3) R*8 O FGS fine lock vector (unit vector) c 1st index: vector component c C FGS_IMA_VEC(3) R*8 O Distortion corrected image space c star vector: c C FGS_OBJ_VEC(3) R*8 O Distortion corrected object space c star vector c C FGS_VEH_VEC(3) R*8 O Star Vector in ST vehicle space C Some local variables: C --------------------- C THETA R*8 Angular component of direction c cosine OF Star vector (radians) C R R*8 Radial component of direction c cosine of Star Vector (radians) C RHOPRIME R*8 Polar radius (radians) to FGS c instantaneous line of sight (LOS) C PHI R*8 Polar angle (radians) to c instantaneous LOS C RPRIME R*8 Polar radius (radians) of the fine c error signal correction term C PSI R*8 Polar angle (radians) of the fine c error signal correction term C DELTA R*8 Spherical excess (radians) C C TERMINOLOGY: C ___________ C Note: terminology and wording below need to be tightened. Also c need to cross-reference with PE, PASS, Flight software, and Astrometry c Team terminology. c FGS (or Star Selector) line of sight: the direction where the c center of the Star Selector is pointing c Guide Star position : position of the star's center as seen in the c FGS. In practice, this is the line of sight defined above, c corrected by the "fine error signal" read from the PMT's. c Guide Star vector: unit vector for the star position c defined above. C DEVELOPMENT HISTORY: C ------------------- C Note: C Original algorithms are from Perkin-Elmer and were implemented c by CSC in PASS. The present code is essentially the PASS code c with: c - minor changes in variable names, c - concatenation of several subroutines, c - additional explanatory comments, c - and inclusion of interfaces to OMS. c The development history below reflects only those changes since c incorporation to OMS. Refer to PASS code for previous history. c c List of PASS modules: c STABERR c UUTATMATD c IM2OBS c FVECTOR c FERRSIG c ROTVEC c OPTDIS c UUTMTMPY8 c UUTQTOMT8 C AUTHOR VERSION DATE CHANGE C ------------ ------- -------- ------------------------ c Pierre Bely 1.0 Sept 87 Original code (see note above) c Lattanzi and c Bucciarelli 2.0 Nov 91 Computation of PHI (see below) c Input data units c c c c============================================================================== IMPLICIT NONE INTEGER FGS_BAY INTEGER*4 PMTXA, PMTXB, PMTYA, PMTYB REAL*8 K1X, K1Y, K0X, K0Y, RX, RY REAL*8 ALEVER, BLEVER, AOFFST, BOFFST, THETA_A, THETA_B INTEGER*4 NX, NY, MX, MY REAL*8 ALPHA(6,6,3), BETA(6,6,3) REAL*8 GM(3) REAL*8 FGS_POS(3),FGS_IMA_VEC(3),FGS_OBJ_VEC(3),FGS_VEH_VEC(3) real*8 pi REAL*8 QX, QY, min2rad, sec2rad, deg2rad REAL*8 X, y, EPS, A, B, C, D, COSA, COSB, # SINA, SINB, THA, THB, SR, CR, SRHO, CRHO, GNU, # THETA, R, SINR,RHOPRIME, RPRIME, PSI, DELTA, # PHI REAL*8 CF, DF(3), XK1, YL1 INTEGER*4 I,J,K, L REAL*8 RTEMP REAL*8 TVS(3,3,3), TVS_MAT(3,3) DATA EPS /13.D-18/ DATA DF /3*-999.0D0/ c------------------------------ Hardcoded for now ----------------- c Data specific to FGS (to be drawn later from PDB's file SCHF) ! ! WRITE(5,*)'FGS_BAY IN SUBROUTINE=',FGS_BAY ! C ALEVER = 411.7494D0 ! C BLEVER = 411.7494D0 ! C AOFFST = 0.0D0 THESE ARE NOW ARGUMENTS CRC 22MAR1994 ! C BOFFST = 0.0D0 ! K0X = 0.0*80.0D-6 ! in arc minutes ! K0Y = 0.0*80D-6 ! C K1X = 256 * 80.0 D-6 ! C K1Y = 256 * 80.0 D-6 ! ! DATA NX,NY,MX,MY/6,6,6,6/ ! c--------------------------------------------------------------------- c conversion factor to radians: pi= dacos(-1.d0) sec2rad = DACOS(-1.D0) / (180.D0*60.D0*60.D0) min2rad = DACOS(-1.D0) / (180.D0*60.D0) deg2rad = DACOS(-1.D0)/180.D0 c-------------------------------------------------------------- C COMPUTE FGS FINE ERROR SIGNAL c-------------------------------------------------------------- c To avoid influence of flux variations, determine the star c position in the detector space not by the straight A-B but c normalizing to the instantaneous total flux: c error signal = (A-B)/(A+B) c Compute normalized x and y error signal (RX and RY) c avoid cases with no counts: if (PMTXA+PMTXB .gt. 0.001D0) then QX=DBLE((PMTXA-PMTXB))/DBLE((PMTXA+PMTXB)) RX = (K0X+K1X*QX)*min2rad else RX=0. endif if (PMTYA+PMTYB .gt. 0.001D0) then QY =DBLE((PMTYA-PMTYB))/DBLE((PMTYA+PMTYB)) RY = (K0Y+K1Y*QY)*min2rad else RY=0. endif C c---------------------------------------------------------------------- c COMPUTE FINE-LOCK FGS VECTOR USING ERROR SIGNALS c---------------------------------------------------------------------- c For convenience find sine and cosine of Star Selector lever arms A and B: COSA = DCOS(ALEVER*min2rad) COSB = DCOS(BLEVER*min2rad) SINA = DSIN(ALEVER*min2rad) SINB = DSIN(BLEVER*min2rad) c Correct the Star Selector angles for zero offset, and convert to radians THA = (THETA_A + AOFFST)*sec2rad THB = (THETA_B + BOFFST)*sec2rad C 1. Determine the polar coordinates of the FGS (star selector) c line of sight (rho' and Phi in Perkin Elmer notation) c a. First determine the radius (rho'): c compute a value for local variable A and dont let A c have an absolute value greater than 1.0 c (note that A is just a calculation convenience and has c nothing to do with the Star Selector arm A) A = COSA*COSB-SINA*SINB*DCOS(THB-THA) IF(DABS(A).GE.1.D0+EPS) A=DSIGN(1.D0,A) RHOPRIME = DACOS(A) SRHO = DSIN(RHOPRIME ) CRHO = DCOS(RHOPRIME ) c b. Then, determine argument (phi): c Note that value X is only for the calculation and has nothing to do c with coordinates. c compute a value for local variable X and dont let X c have an absolute value greater than 1.0 c ***** NOTE: the original PASS computation has been changed ******* c ***** (Nov 18, 1991) c X = (COSB-COSA*CRHO)/(SINA*SRHO) c c c IF(DABS(X).GE.1.D0+EPS) X=DSIGN(1.D0,X) c c c PHI = THA+DACOS(X) c X = -sin(tha-thb) Y = cosa*cos(tha-thb)+sina*cosb/sinb ! Replace with standard DATAN2 function PHI = tha+arctan(x,y,pi) PHI = tha + DATAN2(x,y) c 2. Determine the polar coordinates (r' and Psi in P-E c nomenclature) of the star position in the c detector frame (centered on the Star Selector line of sight) RPRIME = SQRT(RX **2+RY **2) PSI = 0.D0 IF((RY .NE.0.D0).AND.(RX .NE.0.D0)) PSI = ATAN2(RY,RX) c For convenience calculate also the sine and cosine of r' CR = DCOS(RPRIME) SR = DSIN(RPRIME) c 3. Then, find the polar coordinates of the star position in the FGS c image space. c Note that this calculation must be done in the FGS image space c not in the object space - otherwise the spherical excess c correction would be erroneous (see PE document PR-841 - Appendix B) c c a. First, calculate the spherical excess (delta) of the Star c Selector triangle: c Note that B is just a calculation convenience and has nothing to c do with the star selector lever arm B c compute a value for local variable B and dont let B c have an absolute value greater than 1.0 B = DCOS(PHI -THA)*DCOS(THB-THA)+ & DSIN(PHI -THA)*DSIN(THB-THA)*COSA IF(DABS(B).GE.1.D0+EPS) B=DSIGN(1.D0,B) c nu (GNU) is the K angle in the PJK star selector analysis triangle GNU = DACOS(B) DELTA = PHI -THB+GNU c b. Then find the radius to star (R) - (in image space) c note that C is just a calculation convenience c Compute a value for local variable C and dont let C c have an absolute value greater than 1.0 C = CRHO*CR-SRHO*SR*DCOS(PSI +DELTA -PHI ) IF(DABS(C).GE.1.D0+EPS) C=DSIGN(1.D0,C) R = DACOS(C) c c. Finally calculate the argument: c Again, D is just a convenience variable: c Compute a value for local variable D and dont let D c have an absolute value greater than 1.0 D = (SR*DSIN(PSI +DELTA -PHI ))/DSIN(R ) IF(DABS(D).GE.1.D0+EPS) D=DSIGN(1.D0,D) THETA = PHI +DASIN(D) c 3. Derive the components of the Star vector in FGS image c space (optically distorted): c Define SINR as SINR = DSIN(R) c Noting that the projection of the R arc onto the c xy plane is sinR, the vector components X and Y are: FGS_POS(1) = SINR*DCOS(THETA) FGS_POS(2) = SINR*DSIN(THETA) c The Z component is determined from (unit vector): FGS_POS(3) = DSQRT(1.-(FGS_POS(1)**2+FGS_POS(2)**2)) c c--------------------------------------------------------------- c CORRECT FOR OPTICAL DISTORTION c---------------------------------------------------------------- c c Note that this correction is to be done using the c "distorted" to "true" transform polynomial coefficients ! c Note also that it is necessary to go from c - polar coordinates in distorted image space c - to cartesian coordinates in image space, c since the correction transformation is in cartesian coordinates c (It may have been more logical to have it polar coordinates ?) c but then to go back to c - polar coordinates in optically corrected image space c to apply the optical magnification (conserves the argument) c (done in following section), and finally to c - cartesian coordinates in object space c to obtain the direction cosines. C Compute the distortion correction factor for the x-component CF=0.0D0 DO K=1,NX DO L=1,NY c avoid 0**0 uncertainty if((FGS_POS(1).eq.0.0d0) .and. (k.eq.1)) then XK1 = 1.0D0 else XK1 = FGS_POS(1)**(K-1) endif if((FGS_POS(2).eq.0.0d0) .and. (l.eq.1)) then YL1 = 1.0D0 else YL1 = FGS_POS(2)**(L-1) endif CF=CF+ALPHA(K,L,FGS_BAY)*XK1*YL1 enddo enddo cc CC subtract the correction factor from the x-component ! write(15,*)'FGS_POS1',FGS_POS(1),CF FGS_IMA_VEC(1)=FGS_POS(1)-CF CC CC compute the distortion correction factor for the y-component CF=0.0D0 DO K=1,MX DO L=1,MY IF((FGS_POS(1).EQ.0.0D0) .AND. (K.EQ.1)) THEN XK1 = 1.0D0 ELSE XK1 = FGS_POS(1)**(K-1) ENDIF IF((FGS_POS(2).EQ.0.0D0) .AND. (L.EQ.1)) THEN YL1 = 1.0D0 ELSE YL1 = FGS_POS(2)**(L-1) ENDIF CF=CF+BETA(K,L,FGS_BAY)*XK1*YL1 ENDDO ENDDO CC CC subtract the correction factor from the y-component ! write(15,*)'FGS_POS2',FGS_POS(1),CF FGS_IMA_VEC(2)=FGS_POS(2)-CF CC CC calculate the distortion corrected z-component FGS_IMA_VEC(3)= & DSQRT(1.0D0-(FGS_IMA_VEC(1)**2)-(FGS_IMA_VEC(2)**2)) c--------------------------------------------------------------------- c CONVERT FROM IMAGE SPACE TO OBJECT SPACE c--------------------------------------------------------------------- C c See remarks above: go to polar coordinates to apply c the magnification (conserves the argument), c then return to cartesian RTEMP=DSQRT((FGS_IMA_VEC(1)**2)+(FGS_IMA_VEC(2)**2)) c Object to image ratio is sin(R)/sin(R') with c Rtemp=sin(R') and R=R'/magnification c Ref PE report PR 841, Figure 2-12 FGS_OBJ_VEC(1)=(FGS_IMA_VEC(1)/RTEMP)* & DSIN((DASIN(RTEMP))/GM(FGS_BAY)) FGS_OBJ_VEC(2)=(FGS_IMA_VEC(2)/RTEMP)* & DSIN((DASIN(RTEMP))/GM(FGS_BAY)) FGS_OBJ_VEC(3)=DSQRT(1.0D0-(FGS_OBJ_VEC(1)**2)- & (FGS_OBJ_VEC(2)**2)) c-------------------------------------------------------------------------- c CONVERT FGS VECTOR TO VEHICLE SPACE c-------------------------------------------------------------------------- c call MATMUL routine to multiply matrices so as to convert c the FGS object space vector to the ST vehicle reference c frame. c ********** This subroutine is unused here ************* c c Let TVS_MAT be the current transformation matrix. It is defined as: do i=1,3 do j=1,3 TVS_MAT(i,j)= TVS(i,j,FGS_bay) enddo enddo ! CALL MATMULT(TVS_MAT, FGS_OBJ_VEC, FGS_VEH_VEC, 3,3,1) !unused FGS_VEH_VEC(1)=TVS_MAT(1,1)*FGS_OBJ_VEC(1)+ & TVS_MAT(1,2)*FGS_OBJ_VEC(2)+TVS_MAT(1,3)*FGS_OBJ_VEC(3) FGS_VEH_VEC(2)=TVS_MAT(2,1)*FGS_OBJ_VEC(1)+ & TVS_MAT(2,2)*FGS_OBJ_VEC(2)+TVS_MAT(2,3)*FGS_OBJ_VEC(3) FGS_VEH_VEC(3)=TVS_MAT(3,1)*FGS_OBJ_VEC(1)+ & TVS_MAT(3,2)*FGS_OBJ_VEC(2)+TVS_MAT(3,3)*FGS_OBJ_VEC(3) FGS_VEH_VEC(1)=FGS_VEH_VEC(1)/sec2rad FGS_VEH_VEC(2)=FGS_VEH_VEC(2)/sec2rad FGS_VEH_VEC(3)=FGS_VEH_VEC(3)/sec2rad c---------------------------------------------------------------------------- return end ******************************************************************************