C @(#)atref.for 17.1.1.1 (ES0-DMD) 01/25/02 17:55:13 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C @(#)atref.for 17.1.1.1 (ESO-IPG) 17:55:13 01/25/02 PROGRAM ATREF C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.IDENTIFICATION: program ATREF C.PURPOSE: Compensation for atmospheric refraction effects by modification C of the nominal "OPTOPUS" starplate (X,Y) coordinate values. C.AUTHOR: A. Gemmo Padova Department of Astronomy C.VERSION: 860414 G. Lund, P. Angebault Version n.10 for HP machine. C 910412 A. Gemmo Version running in MIDAS env. C 910603 A. Gemmo Version for OPTOPUS context. C 910906 A. Gemmo Modified to use TBLSER routine. C------------------------------------------------------------------------------- IMPLICIT NONE C C *** C INTEGER MADRID(1) INTEGER NRINP INTEGER NCINP INTEGER NROUT INTEGER NCOUT INTEGER NCPAR INTEGER NCOL PARAMETER (NCINP=6) PARAMETER (NCOUT=6) PARAMETER (NCPAR=3) PARAMETER (NROUT=300) INTEGER ISTAT,IAC,IAV INTEGER KUN,KNUL INTEGER NACINP,NARINP,NSINP INTEGER OUTTYP,OUTCOL(10),COLNUM(10) INTEGER TIDINP INTEGER TIDOUT INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9,J INTEGER IMONTH,IMON1 INTEGER IO,ITIME INTEGER INTD,ITST,ITEN,ID,IE,ID1,IE1 INTEGER IWAVL INTEGER IEXPTI,IRANGE,ILAM1,ILAM2 INTEGER ISTSTH,ISTSTM,ISTENH,ISTENM INTEGER IUTSTH,IUTSTM,IUTENH,IUTENM INTEGER IBETST,IBETEN INTEGER DUN,STAT C C *** C DOUBLE PRECISION X(NROUT),Y(NROUT),Z(NROUT) DOUBLE PRECISION PI DOUBLE PRECISION CRA,CDEC,CRARAD,CDECRA DOUBLE PRECISION DATE DOUBLE PRECISION ST1(12),ST2(12) DOUBLE PRECISION SSLOT,ESLOT,ST,SST DOUBLE PRECISION AHS,AHE DOUBLE PRECISION TESS,TESE DOUBLE PRECISION ERINT,SERR DOUBLE PRECISION AH(NROUT),COZ(NROUT),ZZ(NROUT),ERR(NROUT) DOUBLE PRECISION ZIN(3),ZOS(3) DOUBLE PRECISION OBLAT,ZIPHI,TAPHI,A(3),BET(3),ABET(3),SIG DOUBLE PRECISION ZIO,AHIO,SS,TIME,TTIME,BETA2 DOUBLE PRECISION EXPTI,EXPTI1,SEXPT,SEX,EXPST,EXPEN DOUBLE PRECISION STSTH,STENH,STSTM,STENM,D,EXSTUT,EXENUT DOUBLE PRECISION UTSTH,UTENH,UTSTM,UTENM DOUBLE PRECISION EXST,EXEN,ZOST,ZOEN,ZIST,ZIEN DOUBLE PRECISION ABETST,ABETEN,TEST,BETST,BETEN,RANGE DOUBLE PRECISION LAMBD1,LAMBD2,LA1(22),LA2(22) DOUBLE PRECISION CHINT,CHINA,FRAC,CHERR,CORRX,CORRY DOUBLE PRECISION DX(NROUT),DY(NROUT),XO(NROUT),YO(NROUT) DOUBLE PRECISION STD,STF,WAVL DOUBLE PRECISION YEAR,MONTH,DAY,EPP C C *** C CHARACTER*80 STRING(20) CHARACTER*60 INPFIL CHARACTER*60 OUTFIL CHARACTER*16 LABEL(NCPAR),OUTLAB CHARACTER*16 UNIT(NCPAR),OUTUNI CHARACTER*16 OUTFOR CHARACTER*16 FORMC16,FORMC1,FORMR8,FORMI4 CHARACTER*16 IDENT(NROUT) CHARACTER*1 TYPE(NROUT) CHARACTER*1 CHECK(NROUT) CHARACTER*1 ASTFLA C C *** C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/ MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C *** C DATA LABEL/'X ','Y ','Z '/ DATA UNIT /'MICRONS ','MICRONS ','MICRONS '/ C C *** C DATA FORMC16/'A16'/ DATA FORMC1 /'A1' / DATA FORMR8 /'F7.0'/ DATA FORMI4 /'I3'/ DATA PI/3.14159265358979D0/ C C *** sidereal time data C DATA ST1/3.46,5.19,6.51,7.93,9.36,11.29, 2 13.32,14.44,17.69,19.83,22.17,25.12/ DATA ST2/10.27,12.92,15.05,17.39,19.63,21.76,24.0, 2 26.14,27.66,28.88,30.1,31.83/ C C *** wavelength data C DATA LA1/3800.,4000.,4200.,4400.,4600.,4800.,5000., 2 5200.,5400.,5600.,5800.,6000.,6200.,6400., 3 6600.,6800.,7000.,7200.,7400.,7600.,7800., 4 8000./ DATA LA2/1.23,1.02,.83,.69,.58,.49,.4,.33,.26,.19,.14, 2 .09,.05,.02,-.02,-.06,-.1,-.125,-.15,-.175,-.195,-.21/ C C *** start the code C CALL STSPRO('ATREF') C C *** get the input table and the output table C CALL STKRDC('INPUTFI',1,1,60,IAC,INPFIL,KUN,KNUL,ISTAT) CALL STKRDC('OUTPUTF',1,1,60,IAC,OUTFIL,KUN,KNUL,ISTAT) C C *** open input table C CALL TBTOPN(INPFIL,F_I_MODE,TIDINP,ISTAT) C C *** get information about input table C CALL TBIGET(TIDINP,NCOL,NRINP,NSINP,NACINP,NARINP,ISTAT) C C *** create the output table C CALL TBTINI(OUTFIL,0,F_O_MODE,NCOUT,NROUT,TIDOUT,ISTAT) C OUTTYP = D_C_FORMAT !initialize ident column OUTFOR = FORMC16 OUTUNI = ' ' OUTLAB = 'IDENT ' CALL TBCINI(TIDOUT,OUTTYP,16,OUTFOR,OUTUNI, 2 OUTLAB,OUTCOL(1),ISTAT) C OUTTYP = D_C_FORMAT !initialize type column OUTFOR = FORMC1 OUTUNI = ' ' OUTLAB = 'TYPE ' CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI, 2 OUTLAB,OUTCOL(2),ISTAT) C DO I1=1,NCPAR OUTTYP = D_R8_FORMAT !initialize x y z columns OUTFOR = FORMR8 OUTUNI = UNIT(I1) OUTLAB = LABEL(I1) CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI, 2 OUTLAB,OUTCOL(2+I1),ISTAT) ENDDO C OUTTYP = D_I4_FORMAT !initialize number column OUTFOR = FORMI4 OUTUNI = ' ' OUTLAB = 'NUMBER ' CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI, 2 OUTLAB,OUTCOL(6),ISTAT) C C *** read center coordinates C CALL STKRDD('PLATEC1',1,1,IAV,CRA,KUN,KNUL,ISTAT) CALL STKRDD('PLATEC1',2,1,IAV,CDEC,KUN,KNUL,ISTAT) C C *** do the work C IF(CDEC.LT.-89.9)THEN WRITE(STRING(2),922) 922 FORMAT('*** FATAL: DEC of center can`t be lower than -89.9 deg 2rees!!') CALL STETER(9,STRING(2)) ELSE IF(CDEC.GT.30.0)THEN WRITE(STRING(2),923) 923 FORMAT('*** FATAL: DEC of center can`t be higher than +30.0 de 2grees!!') CALL STETER(9,STRING(2)) ENDIF C C *** convert CRA and CDEC to radians C CRARAD = CRA*1.5D1*(PI/1.8D2) CDECRA = CDEC*(PI/1.8D2) C C *** read year, month and day of the observation C CALL STKRDD('DATE',1,1,IAV,YEAR,KUN,KNUL,ISTAT) CALL STKRDD('DATE',2,1,IAV,MONTH,KUN,KNUL,ISTAT) CALL STKRDD('DATE',3,1,IAV,DAY,KUN,KNUL,ISTAT) C IF(YEAR.NE.0..AND.MONTH.EQ.0..AND.DAY.EQ.0.)THEN EPP = YEAR CALL DAYMON(EPP,YEAR,MONTH,DAY) ELSE IF(YEAR.NE.0..AND.MONTH.EQ.0..AND.DAY.NE.0.)THEN CALL STETER(9, 2 '*** FATAL: MONTH can`t be 0 if date not in decimals of year!') ELSE IF(YEAR.NE.0..AND.MONTH.NE.0..AND.DAY.EQ.0.)THEN CALL STETER(9, 2 '*** FATAL: DAY can`t be 0 if date not in decimals of year!') ENDIF C C *** determination of hours of darkness (ST) C DATE = MONTH+DAY/30.42 IMONTH = INT(DATE) IF(IMONTH.EQ.12)THEN IMON1 = 1 STD = ST1(IMONTH)+(24.+ST1(IMON1)-ST1(IMONTH))*(DATE-IMONTH) STF = ST2(IMONTH)+(24.+ST2(IMON1)-ST2(IMONTH))*(DATE-IMONTH) ELSE IF(IMONTH.NE.12)THEN IMON1 = IMONTH+1 STD = ST1(IMONTH)+(ST1(IMON1)-ST1(IMONTH))*(DATE-IMONTH) STF = ST2(IMONTH)+(ST2(IMON1)-ST2(IMONTH))*(DATE-IMONTH) ENDIF IF(STD.GT.24.)STD=STD-24. IF(STF.GT.24.)STF=STF-24. WRITE(STRING(3),221) 221 FORMAT(1X,' ') CALL STTPUT(STRING(3),ISTAT) WRITE(STRING(3),222)STD,STF CALL STTPUT(STRING(3),ISTAT) CALL STDWRD(TIDOUT,'NIGHT_BEG',STD,1,1,DUN,STAT) CALL STDWRD(TIDOUT,'NIGHT_END',STF,1,1,DUN,STAT) C C *** read earliest and latest ST at which the exposure could begin and end C CALL STKRDD('STSLOT',1,1,IAV,SSLOT,KUN,KNUL,ISTAT) CALL STKRDD('STSLOT',2,1,IAV,ESLOT,KUN,KNUL,ISTAT) C C *** read exposure time C CALL STKRDD('EXPTIME',1,1,IAV,EXPTI,KUN,KNUL,ISTAT) C SS = ESLOT-SSLOT IF(SS.LT.0.)SS=SS+24. EXPTI1 = EXPTI/60. C IF((EXPTI1).GE.SS)CALL STETER 1 (9,'*** FATAL: Exposure time can`t be > or = st_slot!') C C *** C AHS = 15.*(SSLOT-CRA) IF(AHS.GT.180.)THEN AHS = AHS-360. ELSE IF(AHS.LT.-180.)THEN AHS = 360.+AHS ENDIF C AHE = 15.*(ESLOT-CRA) IF(AHE.GT.180.)THEN AHE = AHE-360. ELSE IF(AHE.LT.-180.)THEN AHE = 360.+AHE ENDIF C TESS = DABS(AHS) TESE = DABS(AHE) IF(TESS.GE.70.)THEN WRITE(STRING(4),944)AHS 944 FORMAT('*** FATAL: START hour angle can`t be >+/-70 deg, and y 2ours is: ',F6.2,' deg') CALL STETER(9,STRING(4)) ENDIF IF(TESE.GE.70.)THEN WRITE(STRING(4),945)AHE 945 FORMAT('*** FATAL: END hour angle can`t be >+/-70 deg, and you 2rs is: ',F6.2,' deg') CALL STETER(9,STRING(4)) ENDIF IF(TESS.GT.60.)THEN WRITE(STRING(4),946)AHS 946 FORMAT('START hour angle is: ',F6.2,' deg, hence approaching t 2elesc. limit = +/-70 deg') CALL STTPUT(STRING(4),ISTAT) ENDIF IF(TESE.GT.60.)THEN WRITE(STRING(4),947)AHE 947 FORMAT('END hour angle is: ',F6.2,' deg, hence approaching tel 2esc. limit = +/-70 deg') CALL STTPUT(STRING(4),ISTAT) ENDIF C C *** C ERINT = 0. SERR = 0. DO I2=1,NROUT AH(I2) = (AHS+(AHE-AHS)*(I2-1)/99.)*(PI/1.8D2) COZ(I2)= -.4886*DSIN(CDECRA)+(.8725*DCOS(CDECRA)*DCOS(AH(I2))) ZZ(I2) = DACOS(COZ(I2)) ERR(I2)= 43.97*(DTAN(ZZ(I2)+.004363)-TAN(ZZ(I2))) ERINT = ERINT+ERR(I2) ENDDO C C *** read ASTFLAG and ST values C CALL STKRDC('ASTFLAG',1,1,1,IAV,ASTFLA,KUN,KNUL,ISTAT) CALL STKRDD('SIDTIME',1,1,IAV,ST,KUN,KNUL,ISTAT) C 222 FORMAT(1X,' Darkness will begin at ST: ',F5.2,' and end at ST: ', 2 F5.2) 510 FORMAT(1X,' ') 553 FORMAT(1X,' Sidereal time for observation: ',F6.2) 554 FORMAT(1X,' ') 555 FORMAT(1X,' Optimal observing time: ',F6.2) 556 FORMAT(1X,' Optimized observational hour angle in degrees: ',F6.2) 557 FORMAT(1X,' Corresponding distance from the zenith: ',F6.2) 558 FORMAT(1X,' Maximum necessary correction (in arcsecs): ',F6.2) 559 FORMAT(1X,' Direction of correction vectors on starplate: ',F7.2) 560 FORMAT(1X,' (i.e. Perpendicular to the projection of the horizon') 561 FORMAT(1X,' on the starplate at the above determined hour angle)') 5551 FORMAT(1X,' ') C 671 FORMAT(1X,' Chosen length of observation: ',I3,' minutes') 672 FORMAT(1X,' Approx. optimal obs. slot (ST): ',I2,'h ', 2 I2,'m to ',I2,'h ',I2,'m') 673 FORMAT(1X,' Corresponding obser. slot (UT): ',I2,'h ', 2 I2,'m to ',I2,'h ',I2,'m') 674 FORMAT(1X,' Corresp. range of corr. vectors: From',I5, 2 ' to',I5,' deg.') 675 FORMAT(1X,' WARNING!!! The needed correction vector varies') 676 FORMAT(1X,' STRONGLY in DIRECTION during your exposure ...') 677 FORMAT(1X,' i.e. by approximately ',I3,' degrees. ') 678 FORMAT(1X,' BE CAREFUL TO RESPECT THE OPTIMAL OBSERVATION TIME') 679 FORMAT(1X,' SLOT GIVEN ABOVE, WHEN OBSERVING AT THE TELESCOPE') C 6661 FORMAT(1X,' ') 6662 FORMAT(1X,' ') 6664 FORMAT(1X,' ') C 777 FORMAT(1X,' Wavelength range for optimisation: ',I4,' to ', 2 I4,' angstroms') 778 FORMAT(1X,' Optimal correction at wavelength: ',I4,' Angstroms' 2) 779 FORMAT(1X,' Chromatic correction needed in X: ',F5.0, 2' microns ') 780 FORMAT(1X,' Chromatic correction needed in Y: ',F5.0, 2' microns ') C CALL STTPUT(STRING(6),ISTAT) IF(ASTFLA.EQ.'N'.OR.ASTFLA.EQ.'n')THEN IF(ST.GE.SSLOT.AND.ST.LE.ESLOT)THEN SST = ST-SSLOT IF(SST.LT.0.)SST=SST+24. ITIME = INT((SST/SS)*99.+1.) IO = ITIME WRITE(STRING(5),510) CALL STTPUT(STRING(5),ISTAT) WRITE(STRING(5),553)ST CALL STTPUT(STRING(5),ISTAT) CALL STDWRD(TIDOUT,'OPT_SID_TIME',ST,1,1,DUN,STAT) ELSE CALL STETER(9,'*** FATAL: Opt_st must be inside st_slot !!') ENDIF ELSE IF(ASTFLA.EQ.'Y'.OR.ASTFLA.EQ.'y')THEN C C C *** Calculation of the "optimal" observing time, for which the integral C *** differential error will be the same for both of the intervals, i.e. C *** from earliest start to "opt. time", and from "opt. time" to latest C *** termination. C C C DO I3=1,NROUT SERR = SERR+(2.*ERR(I3)) IF(SERR.LE.ERINT)THEN GOTO 444 ELSE IF(SERR.GT.ERINT)THEN IO = I3-1 GOTO 445 ENDIF 444 ENDDO C IO = 50 445 TIME = SSLOT+IO*((SS)/99.) TTIME = TIME IF(TIME.GE.24.)TTIME = TIME-24. ST = TTIME SST = TTIME-SSLOT IF(SST.LT.0.)SST=SST+24. C WRITE(STRING(5),554) CALL STTPUT(STRING(5),ISTAT) WRITE(STRING(5),555)TTIME CALL STTPUT(STRING(5),ISTAT) WRITE(STRING(5),5551) CALL STTPUT(STRING(5),ISTAT) C CALL STDWRD(TIDOUT,'OPT_SID_TIME',TTIME,1,1,DUN,STAT) ENDIF C C *** According to "optimised" time (related to 'IO'), determine the C *** corresponding hour angle, zenith distance, maximum differential C *** error (at plate edge), and correction direction - given by the C *** perpendicular (away from the zenith) of the projection of the C *** horizon onto the starplate. C ZIN(1) = DSIN(ZZ(1)) ZIN(2) = DSIN(ZZ(IO)) ZIN(3) = DSIN(ZZ(NROUT)) ZOS(1) = DCOS(ZZ(1)) ZOS(2) = DCOS(ZZ(IO)) ZOS(3) = DCOS(ZZ(NROUT)) C OBLAT = -29.25*PI/1.8D2 ZIPHI = DSIN(OBLAT) TAPHI = DTAN(OBLAT) C A(1) = AH(1) A(2) = AH(IO) A(3) = AH(NROUT) C DO I4=1,3 IF(ZIN(I4).EQ.0.)THEN BET(I4)=0. ELSE IF(ZIN(I4).NE.0.)THEN ABET(I4) = (.4886+DSIN(CDECRA)*ZOS(I4))/ 2 (DCOS(CDECRA)*ZIN(I4)) SIG=1. IF(A(I4).LT.0.)SIG=-1. BET(I4) = SIG*DACOS(ABET(I4)) ENDIF ENDDO C ZIO = ZZ(IO)/(PI/1.8D2) AHIO = AH(IO)/(PI/1.8D2) BETA2 = BET(2)/(PI/1.8D2) C WRITE(STRING(5),556)AHIO CALL STTPUT(STRING(5),ISTAT) CALL STDWRD(TIDOUT,'HOUR_ANGLE',AHIO,1,1,DUN,STAT) WRITE(STRING(5),557)ZIO CALL STTPUT(STRING(5),ISTAT) CALL STDWRD(TIDOUT,'ZENITH_DIST',ZIO,1,1,DUN,STAT) WRITE(STRING(5),558)ERR(IO) CALL STTPUT(STRING(5),ISTAT) CALL STDWRD(TIDOUT,'MAX_CORR',ERR(IO),1,1,DUN,STAT) WRITE(STRING(5),559)BETA2 CALL STTPUT(STRING(5),ISTAT) CALL STDWRD(TIDOUT,'CORR_VECT_DIR',BETA2,1,1,DUN,STAT) WRITE(STRING(5),560) CALL STTPUT(STRING(5),ISTAT) WRITE(STRING(5),561) CALL STTPUT(STRING(5),ISTAT) C C *** determine the approximate optimal observation slot, in ST and UT, C *** according to the expected exposure length at the telescope. C SEXPT = EXPTI/60. SEX = SEXPT*SST/SS EXPST = ST-SEX IF(EXPST.LE.0.)EXPST = EXPST+24. EXPEN = EXPST+SEXPT IF(EXPEN.GT.24.)EXPEN = EXPEN-24. STSTH = INT(EXPST) STENH = INT(EXPEN) STSTM = INT(60.*(EXPST-STSTH)) STENM = INT(60.*(EXPEN-STENH)) D = ((MONTH-1.)*30.42)+DAY-264. IF(D.LT.0.)D=D+365. INTD = INT(D) EXSTUT=EXPST-0.06575*D+4.7156 EXENUT=EXPEN-0.06575*D+4.7156 IF(EXSTUT.LE.0.)EXSTUT=EXSTUT+24. IF(EXENUT.LE.0.)EXENUT=EXENUT+24. UTSTH = INT(EXSTUT) UTENH = INT(EXENUT) UTSTM = INT(60.*(EXSTUT-UTSTH)) UTENM = INT(60.*(EXENUT-UTENH)) C C *** calculation of corresponding range of correction vector angles C EXST = SST-SEX IF(EXST.LE.0.)EXST=EXST+24. EXEN = EXST+SEXPT IF(EXEN.GT.24.)EXEN=EXEN-24. ITST = INT((EXST/SS)*99.)+1. ITEN = INT((EXEN/SS)*99.)+1. ZOST = DCOS(ZZ(ITST)) ZOEN = DCOS(ZZ(ITEN)) ZIST = DSIN(ZZ(ITST)) ZIEN = DSIN(ZZ(ITEN)) ABETST= (.4886+DSIN(CDECRA)*ZOST)/(DCOS(CDECRA)*ZIST) ABETEN= (.4886+DSIN(CDECRA)*ZOEN)/(DCOS(CDECRA)*ZIEN) SIG = 1. TEST = AH(ITST) IF(TEST.LT.0.)SIG=-1. BETST = (SIG/(PI/1.8D2))*DACOS(ABETST) SIG = 1. TEST = AH(ITEN) IF(TEST.LT.0.)SIG=-1. BETEN = (SIG/(PI/1.8D2))*DACOS(ABETEN) RANGE = DABS(BETEN-BETST) C C *** print out the results C WRITE(STRING(6),6661) CALL STTPUT(STRING(6),ISTAT) IEXPTI = EXPTI WRITE(STRING(6),671)IEXPTI CALL STTPUT(STRING(6),ISTAT) CALL STDWRD(TIDOUT,'EXPTIME',EXPTI,1,1,DUN,STAT) ISTSTH = INT(STSTH) ISTSTM = INT(STSTM) ISTENH = INT(STENH) ISTENM = INT(STENM) WRITE(STRING(6),672)ISTSTH,ISTSTM,ISTENH,ISTENM CALL STTPUT(STRING(6),ISTAT) CALL STDWRD(TIDOUT,'STSLOT',EXPST,1,1,DUN,STAT) CALL STDWRD(TIDOUT,'STSLOT',EXPEN,2,1,DUN,STAT) IUTSTH = INT(UTSTH) IUTSTM = INT(UTSTM) IUTENH = INT(UTENH) IUTENM = INT(UTENM) WRITE(STRING(6),673)IUTSTH,IUTSTM,IUTENH,IUTENM CALL STTPUT(STRING(6),ISTAT) CALL STDWRD(TIDOUT,'UTSLOT',EXSTUT,1,1,DUN,STAT) CALL STDWRD(TIDOUT,'UTSLOT',EXENUT,2,1,DUN,STAT) IBETST = INT(BETST) IBETEN = INT(BETEN) WRITE(STRING(6),674)IBETST,IBETEN CALL STDWRD(TIDOUT,'RANGE_CORR_VECT',BETST,1,1,DUN,STAT) CALL STDWRD(TIDOUT,'RANGE_CORR_VECT',BETEN,2,1,DUN,STAT) WRITE(STRING(6),6662) CALL STTPUT(STRING(6),ISTAT) C IF(RANGE.GT.75.)THEN WRITE(STRING(6),675) CALL STTPUT(STRING(6),ISTAT) WRITE(STRING(6),676) CALL STTPUT(STRING(6),ISTAT) IRANGE = INT(RANGE) WRITE(STRING(6),677)IRANGE CALL STTPUT(STRING(6),ISTAT) WRITE(STRING(6),678) CALL STTPUT(STRING(6),ISTAT) WRITE(STRING(6),679) CALL STTPUT(STRING(6),ISTAT) WRITE(STRING(6),6664) CALL STTPUT(STRING(6),ISTAT) ENDIF C C *** calculation of the optimal wavelength for correction of the chromatic C *** dependancy of differential refraction C CALL STKRDD('LAMBDA',1,1,IAV,LAMBD1,KUN,KNUL,ISTAT) CALL STKRDD('LAMBDA',2,1,IAV,LAMBD2,KUN,KNUL,ISTAT) C ILAM1 = INT(LAMBD1) ILAM2 = INT(LAMBD2) WRITE(STRING(7),777)ILAM1,ILAM2 CALL STTPUT(STRING(7),ISTAT) CALL STDWRD(TIDOUT,'LAMBDA',LAMBD1,1,1,DUN,STAT) CALL STDWRD(TIDOUT,'LAMBDA',LAMBD2,2,1,DUN,STAT) IF(LAMBD1.LT.3800.)LAMBD1=3800. IF(LAMBD2.GT.8000.)LAMBD2=8000. ID = INT((LAMBD1-3800.)*21./4200.+1.) IE = INT((LAMBD2-3800.)*21./4200.+1.) IF(ID.LE.0) ID=1 IF(IE.GT.22)IE=22 ID1 = ID+1 IE1 = IE-1 CHINT = 0.5*(LA2(ID)+LA2(IE)) CHINA = LA2(ID) C DO I5=ID1,IE1 CHINT = CHINT+LA2(I5) ENDDO C DO I6=ID1,IE CHINA = CHINA+2.*LA2(I6) TEST = CHINA-LA2(I6) IF(TEST.GE.CHINT)GOTO 1613 ENDDO C 1613 J = I6-1 FRAC = (CHINT+LA2(I6)+LA2(J)-TEST)/(LA2(I6)+LA2(J)) CHERR = LA2(J)+FRAC*(LA2(I6)-LA2(J)) IWAVL = INT(LA1(J)+FRAC*200.) WAVL = DBLE(IWAVL) C C *** calculate the corresponding correction of guidestar coordinates C *** needed to compensate for the shift at the optimised wavelength C CORRX = 140.06*CHERR*DSIN(BET(2))*DTAN(ZZ(IO)) CORRY = 140.06*CHERR*DCOS(BET(2))*DTAN(ZZ(IO)) C WRITE(STRING(7),778)IWAVL CALL STTPUT(STRING(7),ISTAT) CALL STDWRD(TIDOUT,'OPT_WAVEL',WAVL,1,1,DUN,STAT) WRITE(STRING(7),779)CORRX CALL STTPUT(STRING(7),ISTAT) CALL STDWRD(TIDOUT,'CORR_X',CORRX,1,1,DUN,STAT) WRITE(STRING(7),780)CORRY CALL STTPUT(STRING(7),ISTAT) CALL STDWRD(TIDOUT,'CORR_Y',CORRY,1,1,DUN,STAT) C C *** read input table elements C CALL TBLSER(TIDINP,'IDENT',COLNUM(1),ISTAT) CALL TBLSER(TIDINP,'TYPE',COLNUM(2),ISTAT) CALL TBLSER(TIDINP,'CHECK',COLNUM(3),ISTAT) CALL TBLSER(TIDINP,'X',COLNUM(4),ISTAT) CALL TBLSER(TIDINP,'Y',COLNUM(5),ISTAT) CALL TBLSER(TIDINP,'Z',COLNUM(6),ISTAT) C DO I7=1,NRINP CALL TBERDC(TIDINP,I7,COLNUM(1),IDENT(I7),KNUL,ISTAT) CALL TBERDC(TIDINP,I7,COLNUM(2),TYPE(I7),KNUL,ISTAT) CALL TBERDC(TIDINP,I7,COLNUM(3),CHECK(I7),KNUL,ISTAT) CALL TBERDD(TIDINP,I7,COLNUM(4),X(I7),KNUL,ISTAT) CALL TBERDD(TIDINP,I7,COLNUM(5),Y(I7),KNUL,ISTAT) CALL TBERDD(TIDINP,I7,COLNUM(6),Z(I7),KNUL,ISTAT) ENDDO C C *** application of the corrections to the input file: C *** (1) add the differential corrections to all the coordinates; C *** (2) add the chromatic corrections to the guidestar coordinates. C DO I8=1,NRINP DX(I8) = -1.111D-3*ERR(IO)*(DSIN(BET(2))**2)* 2 (X(I8)+Y(I8)/DTAN(BET(2))) DY(I8) = -1.111D-3*ERR(IO)*(DCOS(BET(2))**2)* 2 (Y(I8)+X(I8)*DTAN(BET(2))) XO(I8) = X(I8)+DX(I8) YO(I8) = Y(I8)+DY(I8) C IF(TYPE(I8).NE.'O')THEN DX(I8) = DX(I8)+CORRX DY(I8) = DY(I8)+CORRY XO(I8) = XO(I8)+CORRX YO(I8) = YO(I8)+CORRY ENDIF ENDDO C C *** fill in the output table C DO I9=1,NRINP CALL TBEWRC(TIDOUT,I9,OUTCOL(1),IDENT(I9),ISTAT) CALL TBEWRC(TIDOUT,I9,OUTCOL(2),TYPE(I9),ISTAT) CALL TBEWRD(TIDOUT,I9,OUTCOL(3),XO(I9),ISTAT) CALL TBEWRD(TIDOUT,I9,OUTCOL(4),YO(I9),ISTAT) CALL TBEWRD(TIDOUT,I9,OUTCOL(5),Z(I9),ISTAT) CALL TBEWRI(TIDOUT,I9,OUTCOL(6),I9,ISTAT) ENDDO C C *** initialize row selection flag C CALL TBSINI(TIDOUT,ISTAT) C C *** close the table C CALL TBTCLO(TIDOUT,ISTAT) C C *** over and out C CALL STSEPI END C------------------------------------------------------------------------------- SUBROUTINE DAYMON(EPP,YEAR,MONTH,DAY) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.PURPOSE: Calculate year, month and day from epoch. C.AUTHOR: Alessandra Gemmo Padova Department of Astronomy C.VERSION: 910906 AG Creation C.------------------------------------------------------------------------------ IMPLICIT NONE C C *** C DOUBLE PRECISION YEAR,MONTH,DAY DOUBLE PRECISION DAY1 DOUBLE PRECISION EPP,NDAYS DOUBLE PRECISION RES0,RES1 C C *** C DATA NDAYS/365./ C C *** determine year C YEAR = INT(EPP) C C *** determine month and day C RES0 = DMOD(YEAR,4.0D0) RES1 = DMOD(YEAR,4.0D2) IF(RES0.EQ.0..AND.RES1.NE.0.)NDAYS = 366. C DAY1 = (EPP - YEAR)*NDAYS C IF(DAY1.GT.0..AND.DAY1.LE.31.)THEN MONTH = 1. DAY = DAY1 - 0. ELSE IF(DAY1.GT.31..AND.DAY1.LE.59.)THEN MONTH = 2. DAY = DAY1 - 31. ELSE IF(DAY1.GT.59..AND.DAY1.LE.90.)THEN MONTH = 3. DAY = DAY1 - 59. ELSE IF(DAY1.GT.90..AND.DAY1.LE.120.)THEN MONTH = 4. DAY = DAY1 - 90. ELSE IF(DAY1.GT.120..AND.DAY1.LE.151.)THEN MONTH = 5. DAY = DAY1 - 120. ELSE IF(DAY1.GT.151..AND.DAY1.LE.181.)THEN MONTH = 6. DAY = DAY1 - 151. ELSE IF(DAY1.GT.181..AND.DAY1.LE.212.)THEN MONTH = 7. DAY = DAY1 - 181. ELSE IF(DAY1.GT.212..AND.DAY1.LE.243.)THEN MONTH = 8. DAY = DAY1 - 212. ELSE IF(DAY1.GT.243..AND.DAY1.LE.273.)THEN MONTH = 9. DAY = DAY1 - 243. ELSE IF(DAY1.GT.273..AND.DAY1.LE.304.)THEN MONTH = 10. DAY = DAY1 - 273. ELSE IF(DAY1.GT.304..AND.DAY1.LE.334.)THEN MONTH = 11. DAY = DAY1 - 304. ELSE IF(DAY1.GT.334..AND.DAY1.LE.366.)THEN MONTH = 12. DAY = DAY1 - 334. ENDIF C RETURN END C-------------------------------------------------------------------------------