← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to stdred/optopus/src/atref.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
stdred/optopus/src/atref.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995-2010 European Southern Observatory (ESO)
 
3
C
 
4
C This program is free software; you can redistribute it and/or 
 
5
C modify it under the terms of the GNU General Public License as 
 
6
C published by the Free Software Foundation; either version 2 of 
 
7
C the License, or (at your option) any later version.
 
8
C
 
9
C This program is distributed in the hope that it will be useful,
 
10
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
C GNU General Public License for more details.
 
13
C
 
14
C You should have received a copy of the GNU General Public 
 
15
C License along with this program; if not, write to the Free 
 
16
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
C MA 02139, USA.
 
18
C
 
19
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
C       Internet e-mail: midas@eso.org
 
21
C       Postal address: European Southern Observatory
 
22
C                       Data Management Division 
 
23
C                       Karl-Schwarzschild-Strasse 2
 
24
C                       D 85748 Garching bei Muenchen 
 
25
C                       GERMANY
 
26
C===========================================================================
 
27
C
 
28
      PROGRAM ATREF
 
29
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C.IDENTIFICATION: program ATREF
 
31
C.PURPOSE: Compensation for atmospheric refraction effects by modification
 
32
C          of the nominal "OPTOPUS" starplate (X,Y) coordinate values.
 
33
C.AUTHOR: A. Gemmo       Padova Department of Astronomy
 
34
C.VERSION: 860414 G. Lund, P. Angebault   Version n.10 for HP machine.
 
35
C          910412 A. Gemmo                Version running in MIDAS env.
 
36
C          910603 A. Gemmo                Version for OPTOPUS context.
 
37
C          910906 A. Gemmo                Modified to use TBLSER routine.
 
38
C 
 
39
C 100616        last modif
 
40
C-------------------------------------------------------------------------------
 
41
      IMPLICIT NONE
 
42
C
 
43
C ***
 
44
C
 
45
      INTEGER           MADRID(1)
 
46
      INTEGER           NRINP
 
47
      INTEGER           NCINP
 
48
      INTEGER           NROUT
 
49
      INTEGER           NCOUT
 
50
      INTEGER           NCPAR
 
51
      INTEGER           NCOL
 
52
      PARAMETER         (NCINP=6)
 
53
      PARAMETER         (NCOUT=6)
 
54
      PARAMETER         (NCPAR=3)
 
55
      PARAMETER         (NROUT=300)
 
56
      INTEGER           ISTAT,IAC,IAV
 
57
      INTEGER           KUN,KNUL
 
58
      INTEGER           NACINP,NARINP,NSINP
 
59
      INTEGER           OUTTYP,OUTCOL(10),COLNUM(10)
 
60
      INTEGER           TIDINP
 
61
      INTEGER           TIDOUT
 
62
      INTEGER           I1,I2,I3,I4,I5,I6,I7,I8,I9,J
 
63
      INTEGER           IMONTH,IMON1
 
64
      INTEGER           IO,ITIME
 
65
      INTEGER           INTD,ITST,ITEN,ID,IE,ID1,IE1
 
66
      INTEGER           IWAVL
 
67
      INTEGER           IEXPTI,IRANGE,ILAM1,ILAM2
 
68
      INTEGER           ISTSTH,ISTSTM,ISTENH,ISTENM
 
69
      INTEGER           IUTSTH,IUTSTM,IUTENH,IUTENM
 
70
      INTEGER           IBETST,IBETEN
 
71
      INTEGER           DUN,STAT
 
72
C
 
73
C ***
 
74
C
 
75
      DOUBLE PRECISION          X(NROUT),Y(NROUT),Z(NROUT)
 
76
      DOUBLE PRECISION          PI
 
77
      DOUBLE PRECISION          CRA,CDEC,CRARAD,CDECRA
 
78
      DOUBLE PRECISION          DATE
 
79
      DOUBLE PRECISION          ST1(12),ST2(12)
 
80
      DOUBLE PRECISION          SSLOT,ESLOT,ST,SST 
 
81
      DOUBLE PRECISION          AHS,AHE
 
82
      DOUBLE PRECISION          TESS,TESE
 
83
      DOUBLE PRECISION          ERINT,SERR
 
84
      DOUBLE PRECISION          AH(NROUT),COZ(NROUT),ZZ(NROUT),ERR(NROUT)
 
85
      DOUBLE PRECISION          ZIN(3),ZOS(3)
 
86
      DOUBLE PRECISION          OBLAT,ZIPHI,TAPHI,A(3),BET(3),ABET(3),SIG
 
87
      DOUBLE PRECISION          ZIO,AHIO,SS,TIME,TTIME,BETA2
 
88
      DOUBLE PRECISION          EXPTI,EXPTI1,SEXPT,SEX,EXPST,EXPEN
 
89
      DOUBLE PRECISION          STSTH,STENH,STSTM,STENM,D,EXSTUT,EXENUT
 
90
      DOUBLE PRECISION          UTSTH,UTENH,UTSTM,UTENM
 
91
      DOUBLE PRECISION          EXST,EXEN,ZOST,ZOEN,ZIST,ZIEN
 
92
      DOUBLE PRECISION          ABETST,ABETEN,TEST,BETST,BETEN,RANGE
 
93
      DOUBLE PRECISION          LAMBD1,LAMBD2,LA1(22),LA2(22)
 
94
      DOUBLE PRECISION          CHINT,CHINA,FRAC,CHERR,CORRX,CORRY
 
95
      DOUBLE PRECISION          DX(NROUT),DY(NROUT),XO(NROUT),YO(NROUT)
 
96
      DOUBLE PRECISION          STD,STF,WAVL
 
97
      DOUBLE PRECISION          YEAR,MONTH,DAY,EPP
 
98
C
 
99
C ***
 
100
C
 
101
      CHARACTER*80      STRING(20)
 
102
      CHARACTER*60      INPFIL
 
103
      CHARACTER*60      OUTFIL
 
104
      CHARACTER*16      LABEL(NCPAR),OUTLAB
 
105
      CHARACTER*16      UNIT(NCPAR),OUTUNI
 
106
      CHARACTER*16      OUTFOR
 
107
      CHARACTER*16      FORMC16,FORMC1,FORMR8,FORMI4
 
108
      CHARACTER*16      IDENT(NROUT)
 
109
      CHARACTER*1       TYPE(NROUT)
 
110
      CHARACTER*1       CHECK(NROUT)
 
111
      CHARACTER*1       ASTFLA
 
112
      CHARACTER      CTEST*16
 
113
C
 
114
C ***
 
115
C      
 
116
      INCLUDE           'MID_INCLUDE:ST_DEF.INC'
 
117
      COMMON            /VMR/ MADRID
 
118
      INCLUDE           'MID_INCLUDE:ST_DAT.INC'
 
119
C
 
120
C ***
 
121
C
 
122
      DATA              LABEL/'X       ','Y       ','Z       '/
 
123
      DATA              UNIT /'MICRONS ','MICRONS ','MICRONS '/
 
124
C
 
125
C ***
 
126
C
 
127
      DATA              FORMC16/'A16'/
 
128
      DATA              FORMC1 /'A1' /
 
129
      DATA              FORMR8 /'F7.0'/
 
130
      DATA              FORMI4 /'I3'/
 
131
      DATA              PI/3.14159265358979D0/
 
132
C
 
133
C *** sidereal time data
 
134
C
 
135
      DATA ST1/3.46,5.19,6.51,7.93,9.36,11.29,
 
136
     2          13.32,14.44,17.69,19.83,22.17,25.12/
 
137
      DATA ST2/10.27,12.92,15.05,17.39,19.63,21.76,24.0,
 
138
     2          26.14,27.66,28.88,30.1,31.83/
 
139
C
 
140
C *** wavelength data
 
141
C
 
142
      DATA LA1/3800.,4000.,4200.,4400.,4600.,4800.,5000.,
 
143
     2          5200.,5400.,5600.,5800.,6000.,6200.,6400.,
 
144
     3          6600.,6800.,7000.,7200.,7400.,7600.,7800.,
 
145
     4          8000./
 
146
      DATA LA2/1.23,1.02,.83,.69,.58,.49,.4,.33,.26,.19,.14,
 
147
     2          .09,.05,.02,-.02,-.06,-.1,-.125,-.15,-.175,-.195,-.21/
 
148
C
 
149
C *** start the code
 
150
C
 
151
      CALL STSPRO('ATREF')
 
152
C
 
153
C *** get the input table and the output table
 
154
C
 
155
      CALL STKRDC('INPUTFI',1,1,60,IAC,INPFIL,KUN,KNUL,ISTAT)
 
156
      CALL STKRDC('OUTPUTF',1,1,60,IAC,OUTFIL,KUN,KNUL,ISTAT)
 
157
C
 
158
C *** open input table
 
159
C
 
160
      CALL TBTOPN(INPFIL,F_I_MODE,TIDINP,ISTAT)
 
161
C
 
162
C *** get information about input table
 
163
C
 
164
      CALL TBIGET(TIDINP,NCOL,NRINP,NSINP,NACINP,NARINP,ISTAT)
 
165
C
 
166
C *** create the output table
 
167
C
 
168
      CALL TBTINI(OUTFIL,0,F_O_MODE,NCOUT,NROUT,TIDOUT,ISTAT)
 
169
C
 
170
      OUTTYP = D_C_FORMAT     !initialize ident column
 
171
      OUTFOR = FORMC16
 
172
      OUTUNI = '        '
 
173
      OUTLAB = 'IDENT   '
 
174
      CALL TBCINI(TIDOUT,OUTTYP,16,OUTFOR,OUTUNI,
 
175
     2                     OUTLAB,OUTCOL(1),ISTAT)
 
176
C
 
177
      OUTTYP = D_C_FORMAT     !initialize type column
 
178
      OUTFOR = FORMC1
 
179
      OUTUNI = '        '
 
180
      OUTLAB = 'TYPE    '
 
181
      CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI,
 
182
     2                     OUTLAB,OUTCOL(2),ISTAT)
 
183
C
 
184
      DO I1=1,NCPAR
 
185
         OUTTYP = D_R8_FORMAT  !initialize x y z columns
 
186
         OUTFOR = FORMR8
 
187
         OUTUNI = UNIT(I1)
 
188
         OUTLAB = LABEL(I1)
 
189
         CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI,
 
190
     2                        OUTLAB,OUTCOL(2+I1),ISTAT)
 
191
        ENDDO
 
192
C
 
193
      OUTTYP = D_I4_FORMAT     !initialize number column
 
194
      OUTFOR = FORMI4
 
195
      OUTUNI = '        '
 
196
      OUTLAB = 'NUMBER  '
 
197
      CALL TBCINI(TIDOUT,OUTTYP,1,OUTFOR,OUTUNI,
 
198
     2                     OUTLAB,OUTCOL(6),ISTAT)
 
199
C
 
200
C *** read center coordinates
 
201
C
 
202
      CALL STKRDD('PLATEC1',1,1,IAV,CRA,KUN,KNUL,ISTAT)
 
203
      CALL STKRDD('PLATEC1',2,1,IAV,CDEC,KUN,KNUL,ISTAT)
 
204
C
 
205
C *** do the work
 
206
C
 
207
      IF(CDEC.LT.-89.9)THEN
 
208
        WRITE(STRING(2),922)
 
209
922     FORMAT('*** FATAL: DEC of center can`t be lower than -89.9 deg
 
210
     2rees!!')
 
211
        CALL STETER(9,STRING(2))
 
212
      ELSE IF(CDEC.GT.30.0)THEN
 
213
        WRITE(STRING(2),923)
 
214
923     FORMAT('*** FATAL: DEC of center can`t be higher than +30.0 de
 
215
     2grees!!')
 
216
        CALL STETER(9,STRING(2))
 
217
      ENDIF
 
218
C
 
219
C *** convert CRA and CDEC to radians
 
220
C
 
221
      CRARAD = CRA*1.5D1*(PI/1.8D2)
 
222
      CDECRA = CDEC*(PI/1.8D2)
 
223
C
 
224
C *** read year, month and day of the observation
 
225
C
 
226
      CALL STKRDD('DATE',1,1,IAV,YEAR,KUN,KNUL,ISTAT)
 
227
      CALL STKRDD('DATE',2,1,IAV,MONTH,KUN,KNUL,ISTAT)
 
228
      CALL STKRDD('DATE',3,1,IAV,DAY,KUN,KNUL,ISTAT)
 
229
C
 
230
      IF(YEAR.NE.0..AND.MONTH.EQ.0..AND.DAY.EQ.0.)THEN
 
231
       EPP = YEAR
 
232
       CALL DAYMON(EPP,YEAR,MONTH,DAY)
 
233
      ELSE IF(YEAR.NE.0..AND.MONTH.EQ.0..AND.DAY.NE.0.)THEN
 
234
       CALL STETER(9,
 
235
     2  '*** FATAL: MONTH can`t be 0 if date not in decimals of year!')
 
236
      ELSE IF(YEAR.NE.0..AND.MONTH.NE.0..AND.DAY.EQ.0.)THEN
 
237
       CALL STETER(9,
 
238
     2  '*** FATAL: DAY can`t be 0 if date not in decimals of year!')
 
239
      ENDIF
 
240
C
 
241
C *** determination of hours of darkness (ST)
 
242
C
 
243
      DATE   = MONTH+DAY/30.42
 
244
      IMONTH = INT(DATE)
 
245
      IF(IMONTH.EQ.12)THEN
 
246
        IMON1 = 1
 
247
        STD   = ST1(IMONTH)+(24.+ST1(IMON1)-ST1(IMONTH))*(DATE-IMONTH)
 
248
        STF   = ST2(IMONTH)+(24.+ST2(IMON1)-ST2(IMONTH))*(DATE-IMONTH)
 
249
      ELSE IF(IMONTH.NE.12)THEN
 
250
        IMON1 = IMONTH+1
 
251
        STD   = ST1(IMONTH)+(ST1(IMON1)-ST1(IMONTH))*(DATE-IMONTH)
 
252
        STF   = ST2(IMONTH)+(ST2(IMON1)-ST2(IMONTH))*(DATE-IMONTH)
 
253
      ENDIF
 
254
      IF(STD.GT.24.)STD=STD-24.
 
255
      IF(STF.GT.24.)STF=STF-24.
 
256
      WRITE(STRING(3),221)
 
257
221   FORMAT(1X,'   ')
 
258
      CALL STTPUT(STRING(3),ISTAT)
 
259
      WRITE(STRING(3),222)STD,STF
 
260
      CALL STTPUT(STRING(3),ISTAT)
 
261
      CALL STDWRD(TIDOUT,'NIGHT_BEG',STD,1,1,DUN,STAT)
 
262
      CALL STDWRD(TIDOUT,'NIGHT_END',STF,1,1,DUN,STAT)
 
263
C
 
264
C *** read earliest and latest ST at which the exposure could begin and end
 
265
C
 
266
      CALL STKRDD('STSLOT',1,1,IAV,SSLOT,KUN,KNUL,ISTAT)
 
267
      CALL STKRDD('STSLOT',2,1,IAV,ESLOT,KUN,KNUL,ISTAT)
 
268
C
 
269
C *** read exposure time
 
270
C
 
271
      CALL STKRDD('EXPTIME',1,1,IAV,EXPTI,KUN,KNUL,ISTAT)
 
272
C
 
273
      SS    = ESLOT-SSLOT
 
274
      IF(SS.LT.0.)SS=SS+24.
 
275
      EXPTI1 = EXPTI/60.
 
276
C
 
277
      IF((EXPTI1).GE.SS)CALL STETER
 
278
     1    (9,'*** FATAL: Exposure time can`t be > or = st_slot!')
 
279
C
 
280
C ***
 
281
C
 
282
      AHS = 15.*(SSLOT-CRA)     
 
283
      IF(AHS.GT.180.)THEN
 
284
        AHS = AHS-360.
 
285
      ELSE IF(AHS.LT.-180.)THEN
 
286
        AHS = 360.+AHS
 
287
      ENDIF
 
288
C
 
289
      AHE = 15.*(ESLOT-CRA)
 
290
      IF(AHE.GT.180.)THEN
 
291
        AHE = AHE-360.
 
292
      ELSE IF(AHE.LT.-180.)THEN
 
293
        AHE = 360.+AHE
 
294
      ENDIF
 
295
C
 
296
      TESS = DABS(AHS)
 
297
      TESE = DABS(AHE)
 
298
      IF(TESS.GE.70.)THEN
 
299
        WRITE(STRING(4),944)AHS
 
300
944     FORMAT('*** FATAL: START hour angle can`t be >+/-70 deg, and y   
 
301
     2ours is: ',F6.2,' deg')
 
302
        CALL STETER(9,STRING(4))
 
303
      ENDIF
 
304
      IF(TESE.GE.70.)THEN
 
305
        WRITE(STRING(4),945)AHE
 
306
945     FORMAT('*** FATAL: END hour angle can`t be >+/-70 deg, and you  
 
307
     2rs is: ',F6.2,' deg')
 
308
        CALL STETER(9,STRING(4))
 
309
      ENDIF
 
310
      IF(TESS.GT.60.)THEN
 
311
        WRITE(STRING(4),946)AHS
 
312
946     FORMAT('START hour angle is: ',F6.2,' deg, hence approaching t
 
313
     2elesc. limit = +/-70 deg')
 
314
        CALL STTPUT(STRING(4),ISTAT)
 
315
      ENDIF
 
316
      IF(TESE.GT.60.)THEN
 
317
        WRITE(STRING(4),947)AHE
 
318
947     FORMAT('END hour angle is: ',F6.2,' deg, hence approaching tel
 
319
     2esc. limit = +/-70 deg')
 
320
        CALL STTPUT(STRING(4),ISTAT)
 
321
      ENDIF
 
322
C
 
323
C *** 
 
324
C
 
325
      ERINT = 0.
 
326
      SERR  = 0.
 
327
      DO I2=1,NROUT
 
328
        AH(I2) = (AHS+(AHE-AHS)*(I2-1)/99.)*(PI/1.8D2)
 
329
        COZ(I2)= -.4886*DSIN(CDECRA)+(.8725*DCOS(CDECRA)*DCOS(AH(I2)))
 
330
        ZZ(I2) = DACOS(COZ(I2))
 
331
        ERR(I2)= 43.97*(DTAN(ZZ(I2)+.004363)-TAN(ZZ(I2)))
 
332
        ERINT  = ERINT+ERR(I2)
 
333
      ENDDO
 
334
C
 
335
C *** read ASTFLAG and ST values
 
336
C
 
337
      CALL STKRDC('ASTFLAG',1,1,1,IAV,ASTFLA,KUN,KNUL,ISTAT)
 
338
      CALL STKRDD('SIDTIME',1,1,IAV,ST,KUN,KNUL,ISTAT)
 
339
C
 
340
222   FORMAT(1X,' Darkness will begin at ST: ',F5.2,' and end at ST: ',
 
341
     2       F5.2)
 
342
510   FORMAT(1X,'   ')
 
343
553   FORMAT(1X,' Sidereal time for observation:                 ',F6.2)
 
344
554   FORMAT(1X,'   ')
 
345
555   FORMAT(1X,' Optimal observing time:                        ',F6.2)
 
346
556   FORMAT(1X,' Optimized observational hour angle in degrees: ',F6.2)
 
347
557   FORMAT(1X,' Corresponding distance from the zenith:        ',F6.2)
 
348
558   FORMAT(1X,' Maximum necessary correction (in arcsecs):     ',F6.2)
 
349
559   FORMAT(1X,' Direction of correction vectors on starplate:  ',F7.2)
 
350
560   FORMAT(1X,' (i.e. Perpendicular to the projection of the horizon')
 
351
561   FORMAT(1X,' on the starplate at the above determined hour angle)')
 
352
5551  FORMAT(1X,'   ')
 
353
C
 
354
671   FORMAT(1X,' Chosen length of observation:      ',I3,' minutes')
 
355
672   FORMAT(1X,' Approx. optimal obs. slot (ST):    ',I2,'h ',
 
356
     2       I2,'m to ',I2,'h ',I2,'m')
 
357
673   FORMAT(1X,' Corresponding obser. slot (UT):    ',I2,'h ',
 
358
     2       I2,'m to ',I2,'h ',I2,'m')
 
359
674   FORMAT(1X,' Corresp. range of corr. vectors:   From',I5,
 
360
     2       '  to',I5,' deg.')
 
361
675   FORMAT(1X,' WARNING!!! The needed correction vector varies')
 
362
676   FORMAT(1X,' STRONGLY in DIRECTION during your exposure ...')
 
363
677   FORMAT(1X,' i.e. by approximately ',I3,' degrees.         ')
 
364
678   FORMAT(1X,' BE CAREFUL TO RESPECT THE OPTIMAL OBSERVATION TIME')
 
365
679   FORMAT(1X,' SLOT GIVEN ABOVE, WHEN OBSERVING AT THE TELESCOPE')
 
366
C
 
367
6661  FORMAT(1X,'   ')
 
368
6662  FORMAT(1X,'   ')
 
369
6664  FORMAT(1X,'   ')
 
370
C
 
371
777   FORMAT(1X,' Wavelength range for optimisation: ',I4,' to ',
 
372
     2  I4,' angstroms')
 
373
778   FORMAT(1X,' Optimal correction at wavelength:  ',I4,' Angstroms'
 
374
     2)
 
375
779   FORMAT(1X,' Chromatic correction needed in X:  ',F5.0,
 
376
     2' microns ')
 
377
780   FORMAT(1X,' Chromatic correction needed in Y:  ',F5.0,
 
378
     2' microns ')
 
379
C
 
380
      CALL STTPUT(STRING(6),ISTAT)
 
381
      IF(ASTFLA.EQ.'N'.OR.ASTFLA.EQ.'n')THEN
 
382
       IF(ST.GE.SSLOT.AND.ST.LE.ESLOT)THEN
 
383
        SST   = ST-SSLOT
 
384
        IF(SST.LT.0.)SST=SST+24.
 
385
        ITIME = INT((SST/SS)*99.+1.)  
 
386
        IO    = ITIME
 
387
        WRITE(STRING(5),510)
 
388
        CALL STTPUT(STRING(5),ISTAT)
 
389
        WRITE(STRING(5),553)ST
 
390
        CALL STTPUT(STRING(5),ISTAT)
 
391
        CALL STDWRD(TIDOUT,'OPT_SID_TIME',ST,1,1,DUN,STAT)
 
392
       ELSE
 
393
        CALL STETER(9,'*** FATAL: Opt_st must be inside st_slot !!')
 
394
       ENDIF
 
395
      ELSE IF(ASTFLA.EQ.'Y'.OR.ASTFLA.EQ.'y')THEN
 
396
C
 
397
C
 
398
C *** Calculation of the "optimal" observing time, for which the integral
 
399
C *** differential error will be the same for both of the intervals, i.e.
 
400
C *** from earliest start to "opt. time", and from "opt. time" to latest
 
401
C *** termination.
 
402
C
 
403
C
 
404
C
 
405
      DO I3=1,NROUT
 
406
        SERR   = SERR+(2.*ERR(I3))
 
407
        IF(SERR.LE.ERINT)THEN
 
408
             GOTO 444
 
409
        ELSE IF(SERR.GT.ERINT)THEN
 
410
             IO = I3-1
 
411
           GOTO 445
 
412
        ENDIF
 
413
444   ENDDO
 
414
C
 
415
      IO = 50   
 
416
 
 
417
445   TIME  = SSLOT+IO*((SS)/99.)
 
418
      TTIME = TIME
 
419
      IF(TIME.GE.24.)TTIME = TIME-24.
 
420
      ST = TTIME
 
421
      SST   = TTIME-SSLOT
 
422
      IF(SST.LT.0.)SST=SST+24.
 
423
C
 
424
      WRITE(STRING(5),554)
 
425
      CALL STTPUT(STRING(5),ISTAT)
 
426
      WRITE(STRING(5),555)TTIME
 
427
      CALL STTPUT(STRING(5),ISTAT)
 
428
      WRITE(STRING(5),5551)
 
429
      CALL STTPUT(STRING(5),ISTAT) 
 
430
C
 
431
      CALL STDWRD(TIDOUT,'OPT_SID_TIME',TTIME,1,1,DUN,STAT)
 
432
      ENDIF
 
433
C
 
434
C *** According to "optimised" time (related to 'IO'), determine the 
 
435
C *** corresponding hour angle, zenith distance, maximum differential 
 
436
C *** error (at plate edge), and correction direction - given by the
 
437
C *** perpendicular (away from the zenith) of the projection of the 
 
438
C *** horizon onto the starplate.
 
439
C
 
440
      ZIN(1) = DSIN(ZZ(1))
 
441
      ZIN(2) = DSIN(ZZ(IO))
 
442
      ZIN(3) = DSIN(ZZ(NROUT))
 
443
      ZOS(1) = DCOS(ZZ(1))
 
444
      ZOS(2) = DCOS(ZZ(IO))
 
445
      ZOS(3) = DCOS(ZZ(NROUT))
 
446
C
 
447
      OBLAT = -29.25*PI/1.8D2
 
448
      ZIPHI = DSIN(OBLAT)
 
449
      TAPHI = DTAN(OBLAT)
 
450
C
 
451
      A(1)  = AH(1)
 
452
      A(2)  = AH(IO)
 
453
      A(3)  = AH(NROUT)
 
454
C
 
455
      DO I4=1,3
 
456
        IF(ZIN(I4).EQ.0.)THEN
 
457
            BET(I4)=0.
 
458
        ELSE IF(ZIN(I4).NE.0.)THEN
 
459
            ABET(I4) = (.4886+DSIN(CDECRA)*ZOS(I4))/
 
460
     2                 (DCOS(CDECRA)*ZIN(I4))
 
461
          SIG=1.
 
462
          IF(A(I4).LT.0.)SIG=-1.
 
463
          BET(I4) = SIG*DACOS(ABET(I4))
 
464
        ENDIF
 
465
      ENDDO
 
466
C
 
467
      ZIO   = ZZ(IO)/(PI/1.8D2)
 
468
      AHIO  = AH(IO)/(PI/1.8D2)
 
469
      BETA2 = BET(2)/(PI/1.8D2)
 
470
C
 
471
      WRITE(STRING(5),556)AHIO
 
472
      CALL STTPUT(STRING(5),ISTAT)
 
473
      CALL STDWRD(TIDOUT,'HOUR_ANGLE',AHIO,1,1,DUN,STAT)
 
474
      WRITE(STRING(5),557)ZIO
 
475
      CALL STTPUT(STRING(5),ISTAT)
 
476
      CALL STDWRD(TIDOUT,'ZENITH_DIST',ZIO,1,1,DUN,STAT)
 
477
      WRITE(STRING(5),558)ERR(IO)
 
478
      CALL STTPUT(STRING(5),ISTAT)
 
479
      CALL STDWRD(TIDOUT,'MAX_CORR',ERR(IO),1,1,DUN,STAT)
 
480
      WRITE(STRING(5),559)BETA2
 
481
      CALL STTPUT(STRING(5),ISTAT)
 
482
      CALL STDWRD(TIDOUT,'CORR_VECT_DIR',BETA2,1,1,DUN,STAT)
 
483
      WRITE(STRING(5),560)
 
484
      CALL STTPUT(STRING(5),ISTAT)
 
485
      WRITE(STRING(5),561)
 
486
      CALL STTPUT(STRING(5),ISTAT) 
 
487
C
 
488
C *** determine the approximate optimal observation slot, in ST and UT,
 
489
C *** according to the expected exposure length at the telescope.
 
490
C
 
491
      SEXPT = EXPTI/60.
 
492
      SEX   = SEXPT*SST/SS
 
493
      EXPST = ST-SEX
 
494
      IF(EXPST.LE.0.)EXPST = EXPST+24.
 
495
      EXPEN = EXPST+SEXPT
 
496
      IF(EXPEN.GT.24.)EXPEN = EXPEN-24.
 
497
      STSTH = INT(EXPST)
 
498
      STENH = INT(EXPEN)
 
499
      STSTM = INT(60.*(EXPST-STSTH))
 
500
      STENM = INT(60.*(EXPEN-STENH))
 
501
      D     = ((MONTH-1.)*30.42)+DAY-264.
 
502
      IF(D.LT.0.)D=D+365.
 
503
      INTD  = INT(D)
 
504
      EXSTUT=EXPST-0.06575*D+4.7156
 
505
      EXENUT=EXPEN-0.06575*D+4.7156
 
506
      IF(EXSTUT.LE.0.)EXSTUT=EXSTUT+24.
 
507
      IF(EXENUT.LE.0.)EXENUT=EXENUT+24.
 
508
      UTSTH = INT(EXSTUT)
 
509
      UTENH = INT(EXENUT)
 
510
      UTSTM = INT(60.*(EXSTUT-UTSTH))
 
511
      UTENM = INT(60.*(EXENUT-UTENH))
 
512
C
 
513
C *** calculation of corresponding range of correction vector angles
 
514
C
 
515
      EXST  = SST-SEX
 
516
      IF(EXST.LE.0.)EXST=EXST+24.
 
517
      EXEN  = EXST+SEXPT
 
518
      IF(EXEN.GT.24.)EXEN=EXEN-24.
 
519
      ITST  = INT((EXST/SS)*99.)+1.
 
520
      ITEN  = INT((EXEN/SS)*99.)+1.
 
521
      ZOST  = DCOS(ZZ(ITST))
 
522
      ZOEN  = DCOS(ZZ(ITEN))
 
523
      ZIST  = DSIN(ZZ(ITST))
 
524
      ZIEN  = DSIN(ZZ(ITEN))
 
525
      ABETST= (.4886+DSIN(CDECRA)*ZOST)/(DCOS(CDECRA)*ZIST)
 
526
      ABETEN= (.4886+DSIN(CDECRA)*ZOEN)/(DCOS(CDECRA)*ZIEN)
 
527
      SIG   = 1.
 
528
      TEST  = AH(ITST)
 
529
      IF(TEST.LT.0.)SIG=-1.
 
530
      BETST = (SIG/(PI/1.8D2))*DACOS(ABETST)
 
531
      SIG   = 1.
 
532
      TEST  = AH(ITEN)
 
533
      IF(TEST.LT.0.)SIG=-1.
 
534
      BETEN = (SIG/(PI/1.8D2))*DACOS(ABETEN)
 
535
      RANGE = DABS(BETEN-BETST)
 
536
C
 
537
C *** print out the results
 
538
C
 
539
      WRITE(STRING(6),6661)
 
540
      CALL STTPUT(STRING(6),ISTAT)
 
541
      IEXPTI = EXPTI
 
542
      WRITE(STRING(6),671)IEXPTI
 
543
      CALL STTPUT(STRING(6),ISTAT)
 
544
      CALL STDWRD(TIDOUT,'EXPTIME',EXPTI,1,1,DUN,STAT)
 
545
      ISTSTH = INT(STSTH)
 
546
      ISTSTM = INT(STSTM)
 
547
      ISTENH = INT(STENH)
 
548
      ISTENM = INT(STENM)
 
549
      WRITE(STRING(6),672)ISTSTH,ISTSTM,ISTENH,ISTENM
 
550
      CALL STTPUT(STRING(6),ISTAT)
 
551
      CALL STDWRD(TIDOUT,'STSLOT',EXPST,1,1,DUN,STAT)
 
552
      CALL STDWRD(TIDOUT,'STSLOT',EXPEN,2,1,DUN,STAT)
 
553
      IUTSTH = INT(UTSTH)
 
554
      IUTSTM = INT(UTSTM)
 
555
      IUTENH = INT(UTENH)
 
556
      IUTENM = INT(UTENM)
 
557
      WRITE(STRING(6),673)IUTSTH,IUTSTM,IUTENH,IUTENM
 
558
      CALL STTPUT(STRING(6),ISTAT)
 
559
      CALL STDWRD(TIDOUT,'UTSLOT',EXSTUT,1,1,DUN,STAT)
 
560
      CALL STDWRD(TIDOUT,'UTSLOT',EXENUT,2,1,DUN,STAT)
 
561
      IBETST = INT(BETST)
 
562
      IBETEN = INT(BETEN)
 
563
      WRITE(STRING(6),674)IBETST,IBETEN
 
564
      CALL STDWRD(TIDOUT,'RANGE_CORR_VECT',BETST,1,1,DUN,STAT)
 
565
      CALL STDWRD(TIDOUT,'RANGE_CORR_VECT',BETEN,2,1,DUN,STAT)
 
566
      WRITE(STRING(6),6662)
 
567
      CALL STTPUT(STRING(6),ISTAT)
 
568
C
 
569
      IF(RANGE.GT.75.)THEN
 
570
        WRITE(STRING(6),675)
 
571
        CALL STTPUT(STRING(6),ISTAT)
 
572
        WRITE(STRING(6),676)
 
573
        CALL STTPUT(STRING(6),ISTAT)
 
574
        IRANGE = INT(RANGE)
 
575
        WRITE(STRING(6),677)IRANGE
 
576
        CALL STTPUT(STRING(6),ISTAT)
 
577
        WRITE(STRING(6),678)
 
578
        CALL STTPUT(STRING(6),ISTAT)
 
579
        WRITE(STRING(6),679)
 
580
        CALL STTPUT(STRING(6),ISTAT)
 
581
        WRITE(STRING(6),6664)
 
582
        CALL STTPUT(STRING(6),ISTAT)
 
583
      ENDIF
 
584
C
 
585
C *** calculation of the optimal wavelength for correction of the chromatic
 
586
C *** dependancy of differential refraction
 
587
C
 
588
      CALL STKRDD('LAMBDA',1,1,IAV,LAMBD1,KUN,KNUL,ISTAT)
 
589
      CALL STKRDD('LAMBDA',2,1,IAV,LAMBD2,KUN,KNUL,ISTAT)
 
590
C
 
591
      ILAM1 = INT(LAMBD1)
 
592
      ILAM2 = INT(LAMBD2)
 
593
      WRITE(STRING(7),777)ILAM1,ILAM2
 
594
      CALL STTPUT(STRING(7),ISTAT)
 
595
      CALL STDWRD(TIDOUT,'LAMBDA',LAMBD1,1,1,DUN,STAT)
 
596
      CALL STDWRD(TIDOUT,'LAMBDA',LAMBD2,2,1,DUN,STAT)
 
597
      IF(LAMBD1.LT.3800.)LAMBD1=3800.
 
598
      IF(LAMBD2.GT.8000.)LAMBD2=8000.
 
599
      ID    = INT((LAMBD1-3800.)*21./4200.+1.)
 
600
      IE    = INT((LAMBD2-3800.)*21./4200.+1.)
 
601
      IF(ID.LE.0) ID=1
 
602
      IF(IE.GT.22)IE=22
 
603
      ID1   = ID+1
 
604
      IE1   = IE-1
 
605
      CHINT = 0.5*(LA2(ID)+LA2(IE))
 
606
      CHINA = LA2(ID)
 
607
C       
 
608
      DO I5=ID1,IE1
 
609
         CHINT = CHINT+LA2(I5)
 
610
      ENDDO
 
611
C
 
612
      DO I6=ID1,IE
 
613
         CHINA = CHINA+2.*LA2(I6)
 
614
         TEST  = CHINA-LA2(I6)
 
615
         IF(TEST.GE.CHINT)GOTO 1613
 
616
      ENDDO
 
617
C
 
618
1613      J     = I6-1
 
619
      FRAC  = (CHINT+LA2(I6)+LA2(J)-TEST)/(LA2(I6)+LA2(J))
 
620
      CHERR = LA2(J)+FRAC*(LA2(I6)-LA2(J))
 
621
      IWAVL = INT(LA1(J)+FRAC*200.)
 
622
      WAVL  = DBLE(IWAVL)
 
623
C
 
624
C *** calculate the corresponding correction of guidestar coordinates
 
625
C *** needed to compensate for the shift at the optimised wavelength
 
626
C
 
627
      CORRX = 140.06*CHERR*DSIN(BET(2))*DTAN(ZZ(IO))
 
628
      CORRY = 140.06*CHERR*DCOS(BET(2))*DTAN(ZZ(IO))
 
629
C
 
630
      WRITE(STRING(7),778)IWAVL
 
631
      CALL STTPUT(STRING(7),ISTAT)
 
632
      CALL STDWRD(TIDOUT,'OPT_WAVEL',WAVL,1,1,DUN,STAT)
 
633
      WRITE(STRING(7),779)CORRX
 
634
      CALL STTPUT(STRING(7),ISTAT)
 
635
      CALL STDWRD(TIDOUT,'CORR_X',CORRX,1,1,DUN,STAT)
 
636
      WRITE(STRING(7),780)CORRY
 
637
      CALL STTPUT(STRING(7),ISTAT)
 
638
      CALL STDWRD(TIDOUT,'CORR_Y',CORRY,1,1,DUN,STAT)
 
639
C
 
640
C *** read input table elements
 
641
C
 
642
         CALL TBLSER(TIDINP,'IDENT',COLNUM(1),ISTAT)
 
643
         CALL TBLSER(TIDINP,'TYPE',COLNUM(2),ISTAT)
 
644
         CALL TBLSER(TIDINP,'CHECK',COLNUM(3),ISTAT)
 
645
         CALL TBLSER(TIDINP,'X',COLNUM(4),ISTAT)
 
646
         CALL TBLSER(TIDINP,'Y',COLNUM(5),ISTAT)
 
647
         CALL TBLSER(TIDINP,'Z',COLNUM(6),ISTAT)
 
648
C
 
649
      DO I7=1,NRINP
 
650
         CALL TBERDC(TIDINP,I7,COLNUM(1),CTEST,KNUL,ISTAT)
 
651
         CALL FT_EOS(CTEST,16,IDENT(I7),ISTAT)
 
652
         CALL TBERDC(TIDINP,I7,COLNUM(2),CTEST,KNUL,ISTAT)
 
653
         TYPE(I7) = CTEST(1:1)
 
654
         CALL TBERDC(TIDINP,I7,COLNUM(3),CTEST,KNUL,ISTAT)
 
655
         CHECK(I7) = CTEST(1:1)
 
656
         CALL TBERDD(TIDINP,I7,COLNUM(4),X(I7),KNUL,ISTAT)
 
657
         CALL TBERDD(TIDINP,I7,COLNUM(5),Y(I7),KNUL,ISTAT)
 
658
         CALL TBERDD(TIDINP,I7,COLNUM(6),Z(I7),KNUL,ISTAT)
 
659
       ENDDO
 
660
C
 
661
C *** application of the corrections to the input file:
 
662
C *** (1) add the differential corrections to all the coordinates;
 
663
C *** (2) add the chromatic corrections to the guidestar coordinates.
 
664
C
 
665
      DO I8=1,NRINP
 
666
         DX(I8) = -1.111D-3*ERR(IO)*(DSIN(BET(2))**2)*
 
667
     2                             (X(I8)+Y(I8)/DTAN(BET(2)))
 
668
         DY(I8) = -1.111D-3*ERR(IO)*(DCOS(BET(2))**2)*
 
669
     2                             (Y(I8)+X(I8)*DTAN(BET(2)))
 
670
         XO(I8) = X(I8)+DX(I8)
 
671
         YO(I8) = Y(I8)+DY(I8)
 
672
C
 
673
         IF(TYPE(I8).NE.'O')THEN
 
674
           DX(I8) = DX(I8)+CORRX
 
675
           DY(I8) = DY(I8)+CORRY
 
676
           XO(I8) = XO(I8)+CORRX
 
677
           YO(I8) = YO(I8)+CORRY
 
678
         ENDIF
 
679
      ENDDO   
 
680
C
 
681
C *** fill in the output table
 
682
C
 
683
      DO I9=1,NRINP
 
684
         CALL TBEWRC(TIDOUT,I9,OUTCOL(1),IDENT(I9),ISTAT)
 
685
         CALL TBEWRC(TIDOUT,I9,OUTCOL(2),TYPE(I9),ISTAT)
 
686
         CALL TBEWRD(TIDOUT,I9,OUTCOL(3),XO(I9),ISTAT)
 
687
         CALL TBEWRD(TIDOUT,I9,OUTCOL(4),YO(I9),ISTAT)
 
688
         CALL TBEWRD(TIDOUT,I9,OUTCOL(5),Z(I9),ISTAT)
 
689
         CALL TBEWRI(TIDOUT,I9,OUTCOL(6),I9,ISTAT)
 
690
      ENDDO
 
691
C
 
692
C *** initialize row selection flag
 
693
C
 
694
      CALL TBSINI(TIDOUT,ISTAT)
 
695
C
 
696
C *** close the table
 
697
C
 
698
      CALL TBTCLO(TIDOUT,ISTAT)
 
699
C
 
700
C *** over and out
 
701
C
 
702
      CALL STSEPI
 
703
      END
 
704
C-------------------------------------------------------------------------------
 
705
      SUBROUTINE DAYMON(EPP,YEAR,MONTH,DAY)
 
706
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
707
C.PURPOSE: Calculate year, month and day from epoch.
 
708
C.AUTHOR: Alessandra Gemmo       Padova Department of Astronomy
 
709
C.VERSION: 910906 AG Creation
 
710
C.------------------------------------------------------------------------------
 
711
      IMPLICIT NONE
 
712
C
 
713
C ***
 
714
C
 
715
      DOUBLE PRECISION     YEAR,MONTH,DAY
 
716
      DOUBLE PRECISION     DAY1
 
717
      DOUBLE PRECISION     EPP,NDAYS
 
718
      DOUBLE PRECISION     RES0,RES1
 
719
C
 
720
C ***
 
721
C
 
722
      DATA                 NDAYS/365./
 
723
C
 
724
C *** determine year
 
725
C
 
726
      YEAR = INT(EPP)
 
727
C
 
728
C *** determine month and day
 
729
C
 
730
      RES0 = DMOD(YEAR,4.0D0)
 
731
      RES1 = DMOD(YEAR,4.0D2)
 
732
      IF(RES0.EQ.0..AND.RES1.NE.0.)NDAYS = 366.
 
733
C
 
734
      DAY1 = (EPP - YEAR)*NDAYS
 
735
C
 
736
      IF(DAY1.GT.0..AND.DAY1.LE.31.)THEN
 
737
       MONTH = 1.
 
738
       DAY = DAY1 - 0.
 
739
      ELSE IF(DAY1.GT.31..AND.DAY1.LE.59.)THEN
 
740
       MONTH = 2. 
 
741
       DAY = DAY1 - 31.
 
742
      ELSE IF(DAY1.GT.59..AND.DAY1.LE.90.)THEN
 
743
       MONTH = 3.
 
744
       DAY = DAY1 - 59.
 
745
      ELSE IF(DAY1.GT.90..AND.DAY1.LE.120.)THEN
 
746
       MONTH = 4.
 
747
       DAY = DAY1 - 90.
 
748
      ELSE IF(DAY1.GT.120..AND.DAY1.LE.151.)THEN
 
749
       MONTH = 5.
 
750
       DAY = DAY1 - 120.
 
751
      ELSE IF(DAY1.GT.151..AND.DAY1.LE.181.)THEN
 
752
       MONTH = 6.
 
753
       DAY = DAY1 - 151.
 
754
      ELSE IF(DAY1.GT.181..AND.DAY1.LE.212.)THEN
 
755
       MONTH = 7.
 
756
       DAY = DAY1 - 181.
 
757
      ELSE IF(DAY1.GT.212..AND.DAY1.LE.243.)THEN
 
758
       MONTH = 8.
 
759
       DAY = DAY1 - 212.
 
760
      ELSE IF(DAY1.GT.243..AND.DAY1.LE.273.)THEN
 
761
       MONTH = 9.
 
762
       DAY = DAY1 - 243.
 
763
      ELSE IF(DAY1.GT.273..AND.DAY1.LE.304.)THEN
 
764
       MONTH = 10.
 
765
       DAY = DAY1 - 273.
 
766
      ELSE IF(DAY1.GT.304..AND.DAY1.LE.334.)THEN
 
767
       MONTH = 11.
 
768
       DAY = DAY1 - 304.
 
769
      ELSE IF(DAY1.GT.334..AND.DAY1.LE.366.)THEN
 
770
       MONTH = 12.
 
771
       DAY = DAY1 - 334.
 
772
      ENDIF
 
773
C
 
774
      RETURN
 
775
      END
 
776
C-------------------------------------------------------------------------------