1
C @(#)redsubs2.for 19.1 (ES0-DMD) 02/25/03 13:28:47
2
C===========================================================================
3
C Copyright (C) 1995 European Southern Observatory (ESO)
5
C This program is free software; you can redistribute it and/or
6
C modify it under the terms of the GNU General Public License as
7
C published by the Free Software Foundation; either version 2 of
8
C the License, or (at your option) any later version.
10
C This program is distributed in the hope that it will be useful,
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
C GNU General Public License for more details.
15
C You should have received a copy of the GNU General Public
16
C License along with this program; if not, write to the Free
17
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge,
20
C Corresponding concerning ESO-MIDAS should be addressed as follows:
21
C Internet e-mail: midas@eso.org
22
C Postal address: European Southern Observatory
23
C Data Management Division
24
C Karl-Schwarzschild-Strasse 2
25
C D 85748 Garching bei Muenchen
27
C===========================================================================
29
SUBROUTINE POLYGON(ARRAY,IFLAG,LABEL)
30
C @(#)redsubs2.for 19.1 (ESO-IPG) 02/25/03 13:28:47
31
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32
C.COPYRIGHT (c) European Southern Observatory 1992
35
C.AUTHOR Andrew T. Young
38
C.PURPOSE more subroutines for REDUCE
44
C-----------------------------------------------------------------------------
47
C scans array for items with NSTR=iflag, and does polygon interpolation.
51
REAL ARRAY, XBAR, YBAR, SLOPE, FRACT, YMAX, YMIN
53
INTEGER IFLAG, LEN, LWORD, JBGN, NIGHT, K, J,
54
1 JMIN, JMAX, JEND, K1, K2
56
INCLUDE 'MID_REL_INCL:obs.inc'
59
DIMENSION ARRAY(MXOBS)
61
CHARACTER A, LABEL*(*), CARD*79
75
CARD(:LEN)=LABEL(:LEN)
76
CARD(LEN+1:LEN+6)=' for '
77
CALL JD2DAT(REAL(2400000.5D0+DTZERO(NIGHT)),CARD(LEN+7:))
85
C find data for current night.
86
IF (NITE(J).NE.NIGHT) GO TO 20
87
IF (NSTR(J).EQ.IFLAG)THEN
104
C JBGN, JEND index the first and last data of ANY kind for the night.
105
C JMIN, JMAX index the first and last IFLAG data for the night.
108
C K is count of data for night N.
110
IF (YMAX.GT.YMIN) THEN
113
CALL PLOT(K,XS,YS,'*')
114
CALL TVN(' Time (decimal days) -->')
115
CALL ASKN('OK to interpolate with a polygon?',A)
117
24 CALL TV('Enter C to replace data with a Constant,')
119
1 CALL TVN(' L to smooth with a Linear fit,')
120
CALL ASKN(' or Q to Quit.',A)
122
C be lazy: fit line, but force it horizontal.
123
CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
125
ELSE IF (A.EQ.'L') THEN
127
CALL TV('Not enough data to fit a line safely.')
128
CALL TVN('Try a constant instead.')
131
CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
137
ARRAY(J)=(DJOBS(J)-XBAR)*SLOPE+YBAR
141
ELSE IF (A.EQ.'Y' .OR. A.EQ.'O') THEN
143
ELSE IF (HELP(A)) THEN
144
CALL TV('Reply NO to smooth data.')
149
ELSE IF (YMAX.EQ.YMIN) THEN
150
CALL TV('Only 1 value available for this night!')
157
CALL TV('NO DATA AVAILABLE FOR THIS NIGHT')
169
ELSE IF (J.GT.JMAX) THEN
172
IF (DJOBS(J).GT.XS(K2)) THEN
178
FRACT=(DJOBS(J)-XS(K1))/(XS(K2)-XS(K1))
179
ARRAY(J)=YS(K1)+(YS(K2)-YS(K1))*FRACT
183
41 CARD(66:78)=CARD(53:65)
184
CARD(53:65)=CARD(40:52)
185
CARD(40:52)=CARD(27:39)
186
CARD(27:39)=CARD(14:26)
187
CARD(14:26)=CARD(:13)
188
CARD(:13)='Interpolated '
191
C fudge to get right plot limits....
194
CALL PLOT(-K,XS,YS,' ')
196
CALL PLOT(-1,DJOBS(J),ARRAY(J),'+')
201
CALL PLOT(K,XS,YS,'*')
202
CALL RTNCON(' Time (decimal days) -->',28)
205
C IF (JBGN.GT.NDATA) RETURN
210
SUBROUTINE DARKFIT(NDET,NDUSED,LABEL)
212
C scans SIGNAL for items with NSTR=-1, models dark level vs. DTMP,
213
C and subtracts dark current for detector number NDET.
216
C C A U T I O N : This routine has *** NOT *** been debugged!
221
REAL XLIMS, YLIMS, SIGMIN, SIGMAX, YMAX,
222
1 YMIN, XBAR, YBAR, SLOPE, XBAR1, YBAR1, SLOPE1, TRY, XBAR2,
224
INTEGER NDET, LEN, LWORD, K, J, JMIN,
225
1 JMAX, K3, K1, I, K2, NIX
227
INCLUDE 'MID_REL_INCL:mbands.inc'
228
C PARAMETER (MBANDS=9)
230
PARAMETER (MXCOLR=3*MBANDS)
231
INTEGER NDUSED(MXCOLR)
233
INCLUDE 'MID_REL_INCL:obs.inc'
234
C This defines SIGNAL, NSTR, KBND, etc.
237
INCLUDE 'MID_REL_INCL:kingfit.inc'
238
INCLUDE 'MID_REL_INCL:dlsq.inc'
240
CHARACTER A1, LABEL*16, CARD*79
242
DIMENSION XLIMS(2),YLIMS(2)
258
CARD(:LEN)=LABEL(:LEN)
259
CARD(LEN+1:)=' Log(DARK) vs. Temperature'
263
C find dark data for current detector.
264
IF (NSTR(J).EQ.-1 .AND. KBND(J).EQ.NDET)THEN
268
IF (SIGNAL(J).GT.0.) THEN
271
CARD='Negative dark reading at MJD ='
272
WRITE(CARD(35:),'(F12.5)') DTZERO(NITE(J))+DJOBS(J)
273
CARD(48:)='in channel'
274
WRITE(CARD(55:),'(I3)') NDET
279
SIGMAX=MAX(SIGMAX,SIGNAL(J))
280
SIGMIN=MIN(SIGMIN,SIGNAL(J))
282
C SYMB is night label:
283
IF (NITE(J).LT.10) THEN
284
SYMB(K)=CHAR(ICHAR('0')+NITE(J))
286
SYMB(K)=CHAR(ICHAR('A')+NITE(J)-1)
292
C do this every time, in case only 1 datum.
298
C K is count of data for current detector.
303
IF (SIGMAX.GT.SIGMIN) THEN
306
CALL PLOT(0,78.,22.,'P')
307
C set different-symbol mode:
308
22 CALL PLOT(0,78.,22.,'D')
311
C plot each night with its own symbol.
312
CALL PLOT(K,XS,YS,SYMB)
313
C set same-symbol mode:
314
CALL PLOT(0,78.,22.,'S')
316
CALL TVN(' Temperature -->')
317
CALL ASKN('OK to fit with 2 lines?',A1)
318
ELSE IF (HELP(A1)) THEN
319
CALL TV('Reply YES only if concave upward')
320
CALL TVN('(high at ends, lower in middle)')
324
CALL RTNCON(' Temperature -->', 20)
328
24 CALL TV('Enter C to replace data with a Constant,')
329
CALL TVN(' L to smooth with a Linear fit,')
330
CALL ASKN(' or Q to Quit.',A1)
333
C be lazy: fit line, but force it horizontal.
334
CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
338
ELSE IF (A1.EQ.'L') THEN
340
CALL TV('Not enough data to fit 1 line safely.')
341
CALL TVN('Try a constant instead.')
344
C fit a simple exponential in temperature:
345
CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
352
C Refine one-exponential fit:
361
C use expected ordinates to weight:
362
W(J)=1./EXP(YBAR+X(1,J)*SLOPE)
365
C subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
366
CALL DLSQ(YP2EXP,K,2, 1, 2, NIX, 1, 0, 1.E-3)
371
C now subtract fit from ALL data with this detector:
374
IF ((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR.
375
1 (NSTR(J).EQ.-1 .AND. KBND(J).NE.NDET)) GO TO 28
377
SIGNAL(J)=SIGNAL(J)-EXP((DTMP(J)-XBAR)*SLOPE+YBAR)
378
IF (NSTR(J).EQ.-1) THEN
379
YMAX=MAX(YMAX,SIGNAL(J))
380
YMIN=MIN(YMIN,SIGNAL(J))
387
ELSE IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN
388
C go on to do two-exponential fit at *****.
390
ELSE IF (A1.EQ.'H')THEN
391
CALL TV('If the plot is generally concave upward,')
392
CALL TVN('and the different nights (different symbols)')
393
CALL TVN('appear to follow the same curve, reply YES.')
394
CALL TV('If the plot is straight or concave downward,')
395
CALL TVN('or some nights are offset, reply NO.')
396
CALL TV('Overlapping points are denoted by "$".')
406
ELSE IF (SIGMAX.EQ.SIGMIN) THEN
407
CALL TV('Only 1 dark value available!')
411
IF (NSTR(J).GT.0 .AND. NDUSED(KBND(J)).EQ.NDET)
412
1 SIGNAL(J)=SIGNAL(J)-SIGMIN
417
CALL TV('NO DATA AVAILABLE -- cannot subtract dark')
423
C Here to fit 2 lines (exponentials), and remove their sum.
428
CALL TV('CAUTION: probably too few points to get a fit.')
429
CALL ASKN('Were there any wild points in the plot above?',A1)
432
CALL TV('Then you should try to fit just one line.')
435
ELSE IF (A1.EQ.'N') THEN
437
ELSE IF (A1.EQ.'H') THEN
438
C get help: show plot again.
446
C fit trial line to bottom 1/3 of data.
447
CALL ROBLIN(XS,YS,K3, XBAR1,YBAR1,SLOPE1)
449
C remove lower trial line, and fit the upper line.
454
C move upper third to lower third:
456
TRY=EXP(YS(K2))-EXP(YBAR1+SLOPE1*(XS(K2)-XBAR1))
471
CALL TV('Sorry; can''t fit 2 lines. Try 1 line.')
475
CALL ROBLIN(XS,YS,K1, XBAR2,YBAR2,SLOPE2)
476
IF (SLOPE2.LT.0. .OR. SLOPE2.LT.SLOPE1) THEN
477
C unphysical result; complain.
478
CALL TV('Sorry; attempted fit is nonsense.')
479
CALL RTNCON('Try again, using 1 line.',24)
485
C Now restore saved data to lower K1:
491
C ***** Polish the result. *****
498
XBAR=(XBAR1+XBAR2)*0.5
502
W(I)=1./(EXP(YBAR1+SLOPE1*(XS(I)-XBAR1)) +
503
1 EXP(YBAR2+SLOPE2*(XS(I)-XBAR2)) )
506
PG(1)=EXP(YBAR1+SLOPE1*(XBAR-XBAR1))
508
PG(3)=EXP(YBAR2+SLOPE2*(XBAR-XBAR2))
512
C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
513
CALL DLSQ(YP2EXP, K, 2, 1, 2, NIX, 1, 0, 1.E-3)
515
IF (P(1).LE.0.D0 .OR. P(3).LE.0.D0) THEN
516
CALL TV('Sorry -- unable to fit two lines. Try one.')
524
YBAR1=LOG(P(1))-SLOPE1*(XBAR-XBAR1)
526
YBAR2=LOG(P(3))-SLOPE2*(XBAR-XBAR2)
529
W(J)=1./(EXP(YBAR1+SLOPE1*(XS(I)-XBAR1)) +
530
1 EXP(YBAR2+SLOPE2*(XS(I)-XBAR2)) )
536
C Subtract the interpolated values from ALL data for this det.
539
IF ((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR.
540
1 (NSTR(J).EQ.-1 .AND. KBND(J).NE.NDET)) GO TO 60
541
FIT=EXP(YBAR1+SLOPE1*(DTMP(J)-XBAR1)) +
542
1 EXP(YBAR2+SLOPE2*(DTMP(J)-XBAR2))
543
SIGNAL(J)=SIGNAL(J)-FIT
544
IF (NSTR(J).NE.-1) GO TO 60
545
C use only dark data to set limits for residual plot.
546
YMAX=MAX(YMAX,SIGNAL(J))
547
YMIN=MIN(YMIN,SIGNAL(J))
550
C Now inspect resids. vs. time:
552
61 CARD(:LEN)=LABEL(:LEN)
553
CARD(LEN+1:)=' Dark RESIDUALS vs. Time'
556
C set plot limits....
558
XLIMS(2)=DJOBS(NDATA)
561
CALL PLOT(0,XLIMS,YLIMS,'L')
563
IF (NSTR(J).NE.-1 .OR. KBND(J).NE.NDET) GO TO 80
564
IF (NITE(J).LT.10) THEN
565
A1=CHAR(ICHAR('0')+NITE(J))
567
A1=CHAR(ICHAR('A')+NITE(J)-1)
569
CALL PLOT(-1,DJOBS(J),SIGNAL(J),A1)
572
CALL PLOT(1,0.,0.,'+')
573
CALL RTNCON(' Time (decimal days) -->',28)
581
C Fits Y = P(1)*EXP(P(2)*Z(1)) + P(3)*EXP(P(4)*Z(1))
587
INCLUDE 'MID_REL_INCL:kingfit.inc'
588
INCLUDE 'MID_REL_INCL:dlsq.inc'
590
DOUBLE PRECISION TERM
593
C Start with first term:
596
C KIP is set by DARKFIT.
602
TERM=EXP(MIN(TERM,30.D0))
604
PART(2)=P(1)*TERM*Z(1)
610
C Here for second term:
616
TERM=EXP(MIN(TERM,30.D0))
618
PART(4)=P(3)*TERM*Z(1)
625
SUBROUTINE DARKPOLY(NDET,NDUSED,LABEL)
627
C scans SIGNAL for items with NSTR=-1, and does polygon interpolation
628
C and subtraction of dark current for detector number NDET.
633
REAL XLIMS, YLIMS, SIGMIN, SIGMAX, SIGA, SIGZ, YMAX,
634
1 YMIN, XBAR, YBAR, SLOPE, FRACT
635
INTEGER NDET, LEN, LWORD, JBGN, NIGHT, K, J,
636
1 JMIN, JMAX, JEND, K1, K2, NIX
638
INCLUDE 'MID_REL_INCL:mbands.inc'
639
C PARAMETER (MBANDS=9)
641
PARAMETER (MXCOLR=3*MBANDS)
642
INTEGER NDUSED(MXCOLR)
644
INCLUDE 'MID_REL_INCL:obs.inc'
645
C This defines NITE, SIGNAL, NSTR, KBND, etc.
648
INCLUDE 'MID_REL_INCL:kingfit.inc'
649
INCLUDE 'MID_REL_INCL:dlsq.inc'
651
CHARACTER A1, LABEL*16, CARD*79
653
DIMENSION XLIMS(2),YLIMS(2)
666
DO 100 NIGHT=1,NIGHTS
670
CARD(:LEN)=LABEL(:LEN)
671
CARD(LEN+1:LEN+10)=' DARK on '
672
CALL JD2DAT(REAL(2400000.5D0+DTZERO(NIGHT)),CARD(LEN+11:))
673
CARD(LWORD(CARD)+3:)='= MJD'
674
WRITE (CARD(LWORD(CARD)+1:),'(F7.0)') DTZERO(NIGHT)
682
C find dark data for current night and detector.
683
IF (NITE(J).NE.NIGHT) GO TO 20
684
c KBND holds NDET for DARK data, remember.
685
IF (NSTR(J).EQ.-1 .AND. KBND(J).EQ.NDET)THEN
691
IF (K.EQ.1) SIGA=SIGNAL(J)
694
SIGMAX=MAX(SIGMAX,SIGNAL(J))
695
SIGMIN=MIN(SIGMIN,SIGNAL(J))
704
C JBGN, JEND index the first and last data of ANY kind for the night.
705
C JMIN, JMAX index the first and last DARK data for the night.
706
C SIGA, SIGZ are the first and last dark LEVELS for the night.
707
C SIGMIN, SIGMAX are the lowest and highest dark LEVELS for the night.
713
C K is count of data for night N.
715
IF (SIGMAX.GT.SIGMIN) THEN
717
CALL PLOT(K,XS,YS,'*')
718
CALL TVN(' TIME (decimal days) -->')
719
CALL ASKN('OK to interpolate with a polygon?',A1)
721
24 CALL TV('Enter C to replace data with a Constant,')
723
1 CALL TVN(' L to smooth with a Linear fit,')
724
CALL ASKN(' or Q to Quit.',A1)
726
C be lazy: fit line, but force it horizontal.
727
CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
731
ELSE IF (A1.EQ.'L') THEN
733
CALL TV('Not enough data to fit a line safely.')
734
CALL TVN('Try a constant instead.')
737
CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
739
ELSE IF (HELP(A1)) THEN
740
CALL TV('Reply NO to smooth data.')
758
C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
759
CALL DLSQ(YPLIN,K,2, 1, 2, NIX, 1, 0, 1.E-3)
767
IF((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR.
768
C real star, not observed with current det...
769
1 (NSTR(J).LE.0 .AND. KBND(J).NE.NDET)) GO TO 28
770
C ...or dark, not observed with current det.
771
C OK, subtract interpolated dark line:
772
SIGNAL(J)=SIGNAL(J)-((DJOBS(J)-XBAR)*SLOPE+YBAR)
773
IF (NSTR(J).EQ.-1) THEN
774
C reset limits for residual plot.
775
YMAX=MAX(YMAX,SIGNAL(J))
776
YMIN=MIN(YMIN,SIGNAL(J))
783
CARD(66:78)=CARD(53:65)
784
CARD(53:65)=CARD(40:52)
785
CARD(40:52)=CARD(27:39)
786
CARD(27:39)=CARD(14:26)
787
CARD(14:26)=CARD(:13)
788
CARD(:13)='RESIDUALS of '
790
C set plot limits....
795
CALL PLOT(0,XLIMS,YLIMS,'L')
798
IF (NSTR(J).NE.-1 .OR. KBND(J).NE.NDET) GO TO 29
799
C plot only dark resids., not star data.
800
CALL PLOT(-1,DJOBS(J),SIGNAL(J),'*')
803
CALL PLOT(1,0.,0.,'+')
804
CALL RTNCON(' Time (decimal days) -->',28)
807
ELSE IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN
812
ELSE IF (K.EQ.1) THEN
813
CALL TV('Only 1 datum available for this night!')
817
IF((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR.
818
1 (NSTR(J).LE.0 .AND. KBND(J).NE.NDET)) THEN
821
SIGNAL(J)=SIGNAL(J)-SIGA
827
CALLTV('NO DATA AVAILABLE FOR THIS NIGHT -- no subtraction')
831
C Do polygon interpolation here.
837
IF (NDUSED(KBND(J)).NE.NDET) GO TO 40
838
C skip observations with other detectors.
840
C Before 1st dark (at JMIN); use it.
841
SIGNAL(J)=SIGNAL(J)-SIGA
842
ELSE IF (J.GT.JMAX) THEN
843
C After last dark (at JMAX); use it.
844
SIGNAL(J)=SIGNAL(J)-SIGZ
847
IF (DJOBS(J).GT.XS(K2)) THEN
848
C Step to next interval.
854
FRACT=(DJOBS(J)-XS(K1))/(XS(K2)-XS(K1))
855
SIGNAL(J)=SIGNAL(J)-(YS(K1)+(YS(K2)-YS(K1))*FRACT)
857
C note that we have no plot of residuals at the DARK data points,
858
C as they are negligible for polygon fit.
862
C IF (JBGN.GT.NDATA) RETURN
867
SUBROUTINE ROBLIN (X, Y, N, XBAR, YBAR, SLOPE)
869
C Fits a robust line to N pairs (X,Y), having medians (Xbar, Ybar)
870
C and SLOPE. Method: see Hoaglin, Mosteller & Tukey, Ch.5.
872
C Does NOT mess up input data; but leaves them sorted on X in XDUM,YDUM.
886
INCLUDE 'MID_REL_INCL:obs.inc'
889
COMMON /SCRAT/ XDUM(MXOBS), YDUM(MXOBS), ZDUM(MXOBS)
893
INTEGER I,J,N3,MID,MIDL,MIDR
894
REAL XL,XR,YL,YR,B,B0,B1,B2,RL,RR,DELTA1,DELTAX,DR0,DR1
897
C ***** BEGIN EXECUTION *****
900
CALL TV('Too many data to sort in available slots.')
901
CALL TV('Expand MXOBS and recompile.')
906
IF (X(2).NE.X(1)) THEN
907
SLOPE=(Y(2)-Y(1))/(X(2)-X(1))
918
CALL TV('ROBLIN can''t fit a line to 0 points!')
919
CALL TVN('BUG in program')
930
CALL SORT2(XDUM,YDUM,N)
935
XBAR=(XDUM(MID)+XDUM(MID+1))*0.5
936
YBAR=(ZDUM(MID)+ZDUM(MID+1))*0.5
943
C Pick size of end thirds:
950
IF (MIDL.EQ.N3/2) THEN
952
XL=(XDUM(MIDL)+XDUM(MIDL+1))*0.5
953
XR=(XDUM(MIDR)+XDUM(MIDR-1))*0.5
961
C Now find y medians for left and right groups.
967
IF (MIDL.EQ.N3/2) THEN
969
YL=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5
976
20 ZDUM(I)=YDUM(N+1-I)
979
IF (MIDL.EQ.N3/2) THEN
981
YR=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5
990
C B0 is initial slope.
992
C ***** BEGIN ITERATIONS *****
994
C Put RESIDS. in zdum:
997
30 ZDUM(I)=YDUM(I)-(YBAR+B*(XDUM(I)-XBAR))
1000
IF (MIDL.EQ.N3/2) THEN
1002
RL=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5
1009
40 ZDUM(I)=YDUM(N+1-I)-(YBAR+B*(XDUM(N+1-I)-XBAR))
1012
IF (MIDL.EQ.N3/2) THEN
1014
RR=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5
1024
IF (B1.EQ.B0) GOTO 100
1028
IF (DR1.EQ.0.) GOTO 100
1029
IF (ABS((DR1-DR0)/DR1).LT.1.E-4) GOTO 100
1030
B2=B1-DR1*(B1-B0)/(DR1-DR0)
1032
IF (ABS(B2-B1).LT.1.E-6*ABS(B1)) GO TO 100
1033
C prepare for next round:
1042
C Current value of b is converged. Recalculate all resids.:
1045
110 ZDUM(I)=YDUM(I)-(YBAR+B*(XDUM(I)-XBAR))
1048
IF (MID.EQ.N/2) THEN
1050
YBAR=(ZDUM(MID)+ZDUM(MID+1))*0.5+YBAR
1059
SUBROUTINE SORT1(X,N)
1061
C Sorts X, in place.
1063
C BOOTHROYD'S SHELL-SORT: COMM.ACM 6,445(1963) ALG.201.
1065
C SEE H.LORIN,"SORTING & SORT SYSTEMS" (ADDISON-WESLEY,1975)
1066
C CH.3 & 11 FOR THEORY AND TESTS.
1091
IF(X(I+M).GE.X(I)) GO TO 20
1100
SUBROUTINE SORT2(X,Y,N)
1102
C Sorts (X,Y) pairs on X, in place.
1104
C BOOTHROYD'S SHELL-SORT: COMM.ACM 6,445(1963) ALG.201.
1106
C SEE H.LORIN,"SORTING & SORT SYSTEMS" (ADDISON-WESLEY,1975)
1107
C CH.3 & 11 FOR THEORY AND TESTS.
1132
IF(X(I+M).GE.X(I)) GO TO 20
1146
C YP for use by DLSQ in linear fits.
1149
C fits: Y = P(1) + P(2)*Z(1)
1154
C Includes for dlsq:
1155
INCLUDE 'MID_REL_INCL:kingfit.inc'
1156
INCLUDE 'MID_REL_INCL:dlsq.inc'
1160
YT = P(1) + P(2)*Z(1)
1172
SUBROUTINE SKYSUB(BANDS,KMIN,KMAX,DTYPE,KDIAPHRAGM,EXTIN)
1174
C scans SIGNAL for items with DTYPE='Y', and does modelling,
1175
C interpolation, and subtraction of sky for each passband.
1181
INCLUDE 'MID_REL_INCL:mbands.inc'
1182
C PARAMETER (MBANDS=9)
1184
INCLUDE 'MID_REL_INCL:obs.inc'
1186
COMMON /SCRAT/ XDUM(MXOBS), YDUM(MXOBS), ZDUM(MXOBS)
1191
INTEGER KMIN, KMAX, KDIAPHRAGM
1192
CHARACTER DTYPE(MXOBS)
1196
REAL COSOB, SINOB, RASUN, DESUN, HASUN, SOLONG, ELMOON, BMOON,
1197
1 SINPHI, COSPHI, ELONG, ELROT, UTROT, STUTZ, TWOPI, ST2UT,
1198
2 TNOON1, TNOON2, PI, DEGRAD, ALAT, XLIMS, YLIMS, ALTS, T,
1199
3 YMAX, YMIN, AIRMAX, SALTS, OBLIQ, SIGMIN,
1200
4 SIGMAX, RISE, SET, DAWN, DUSK, UTRAD, TLST, HAMOON, COSDEL,
1201
5 SINDEL, COSHA, XM, YM, ZM, AZMOON, DELAZ, COSEL, XBAR,
1202
6 YBAR, SLOPE, OLDVAR, ELONGMAX, ELONGMIN, ZENITH
1204
INTEGER NB, LEN, LWORD, JBGN, NIGHT, NTWI,
1205
1 K, J, JMIN, JMAX, JEND, MOONUP, NERROR, NIX, NPTS, NFIT
1207
C COMMONS FOR SPHERICAL TRIG.:
1209
COMMON /CSUN/ COSOB,SINOB,RASUN,DESUN,HASUN,SOLONG,ELMOON,BMOON
1210
COMMON /SPHERE/ SINPHI,COSPHI,ELONG,ELROT,UTROT,STUTZ,TWOPI,ST2UT,
1211
1 TNOON1,TNOON2,PI,DEGRAD,ALAT
1213
PARAMETER (MXCOLR=3*MBANDS)
1215
C commons for star catalog:
1216
INCLUDE 'MID_REL_INCL:mstars.inc'
1217
C PARAMETER (MSTARS=1650)
1219
COMMON /SCATA/ STARS(MSTARS)
1221
C Includes for dlsq:
1222
INCLUDE 'MID_REL_INCL:kingfit.inc'
1223
INCLUDE 'MID_REL_INCL:dlsq.inc'
1226
C Commons for 12 Diaphragms:
1228
INTEGER NP2NDIA(12*MBANDS)
1229
INTEGER NDIA2NP(12,MBANDS), NDIAS,NDIA
1231
REAL XNDIA(12,MBANDS)
1232
COMMON /DIACOM/XNDIA,NP2NDIA,NDIA2NP,NDIAS,NDIA
1234
CHARACTER DIANAM(12)*4
1235
COMMON /DIACH/DIANAM
1240
CHARACTER A1, CARD*80, BANDS(MXCOLR)*8, DATSTR*16
1243
INTEGER LORDER(12),L,NN, JDUSK,JDAWN, NPARMS,NREFDIA
1245
LOGICAL TWILIT, SHORT
1248
EXTERNAL YPLIN, YPSKY
1251
DIMENSION XLIMS(2),YLIMS(2)
1252
DIMENSION ALTS(2), SALTS(5)
1254
C true altitude of Moon at rise & set:
1255
DATA ALTS/-.0145,-.0145/
1258
C true altitude of Sun at limit of astronomical twilight:
1259
DATA SALTS/-.309,0.,0.,0.,-.309/
1266
CALL TV(' SKY subtraction:')
1267
CALL TV('You will see graphs of sky brightness, plotted against')
1268
CALL TVN('time, airmass, and lunar elongation, for each band,')
1269
CALL TVN('on each night separately.')
1270
CALL TV('Then you can choose the method of fitting or')
1271
CALL TVN('interpolating sky brightness, for each case.')
1274
CALL PLOT(0,78.,23.,'P')
1276
DO 1000 NB=KMIN,KMAX
1277
LEN=LWORD(BANDS(NB))
1279
DO 1000 NIGHT=1,NIGHTS
1282
2 CARD(:LEN)=BANDS(NB)
1283
CARD(LEN+1:LEN+10)=' SKY on'
1284
CALL JD2DAT(REAL(2400000.5D0+DTZERO(NIGHT)),CARD(LEN+11:))
1298
C find range of sky data for current night and passband.
1299
IF (NITE(J).NE.NIGHT) GO TO 12
1300
IF (KBND(J).NE.NB) GO TO 10
1301
IF (DTYPE(J).EQ.'Y')THEN
1306
YMAX=MAX(YMAX,SIGNAL(J))
1307
YMIN=MIN(YMIN,SIGNAL(J))
1308
AIRMAX=MAX(AIRMAX,AIRM(J))
1310
IF (K.EQ.1) SIGMIN=SIGNAL(J)
1321
C JBGN, JEND index the first and last data of ANY kind for the night.
1322
C JMIN, JMAX index the first and last SKY data for the night & band.
1323
C SIGMIN, SIGMAX are the first and last sky LEVELS for the night.
1325
C K is count of data for night N.
1327
T=(DTZERO(NIGHT)+INT(DJOBS((JBGN+JEND)/2))-51544.5)/36525.
1329
C STUTZR sets sidereal time for later use by UTMOON & UTSUN.
1330
C Set OBLIQ for SUN (called by UTSUN).
1331
OBLIQ=(23.43929-T*(1.30042E-2+T*(1.639E-7-T*5.036E-7)))*DEGRAD
1335
C Adjust YMAX for possible outliers:
1338
C (Note reversal of Y and X.)
1339
YM=(YS((K+1)/2)+YS(K/2+1))*0.5
1343
C ZM is interquartile range.
1344
YMAX=MIN(YMAX,YM+4.*ZM)
1348
C Take shortcut if required:
1354
C set up plots; first, TIME:
1357
XLIMS(1)=DJOBS(JMIN)
1358
XLIMS(2)=DJOBS(JMAX)
1361
CALL PLOT(0,XLIMS,YLIMS,'L')
1363
C plot moonrise & moonset. First, guarantee RISE < SET:
1366
C UTMOON needs T in Julian centuries from J2000.
1367
T=((2400000.5D0-2451545.D0)+DTZERO(NIGHT)+
1368
1 INT(DJOBS((JBGN+JEND)/2)))/36525.
1369
CALL UTMOON(T,ALTS,1,*21)
1370
C UTMOON returns result in UTROT, in /SPHERE/.
1371
RISE=UTROT+INT(DJOBS((JBGN+JEND)/2))
1372
C we have at most 18 intervals between YMIN & YMAX.
1373
CALL PLOT(-1,RISE,(YMIN+14.*YMAX)/15.,'|')
1374
CALL PLOT(-1,RISE,(YMIN*2.+13.*YMAX)/15.,'|')
1375
CALL PLOT(-1,RISE,YMAX,'R')
1376
IF (DJOBS(JMIN).LT.RISE .AND. DJOBS(JMAX).GT.RISE)
1377
1 CARD(LWORD(CARD)+5:)='R = moonrise'
1378
21 CALL UTMOON(T,ALTS,2,*22)
1379
SET=UTROT+INT(DJOBS((JBGN+JEND)/2))
1380
CALL PLOT(-1,SET,(YMIN+14.*YMAX)/15.,'|')
1381
CALL PLOT(-1,SET,(YMIN*2.+13.*YMAX)/15.,'|')
1382
CALL PLOT(-1,SET,YMAX,'S')
1383
IF (DJOBS(JMIN).LT.SET .AND. DJOBS(JMAX).GT.SET)
1384
1 CARD(LWORD(CARD)+5:)='S = moonset'
1387
C RISE & SET now are correct for Moon. Do SUN:
1390
C Here the argument is DAYS from J2000.
1391
CALL UTSUN(T*36525.,SALTS,1,*23)
1392
DUSK=UTROT+INT(DJOBS((JBGN+JEND)/2))
1393
23 CALL UTSUN(T*36525.,SALTS,5,*24)
1394
DAWN=UTROT+INT(DJOBS((JBGN+JEND)/2))
1396
C DUSK & DAWN are now limits of astronomical twilight.
1402
IF (KBND(J).NE.NB) GO TO 26
1403
IF ( (DAWN.LT.DUSK .AND.
1404
1 DJOBS(J).GT.DAWN .AND. DJOBS(J).LT.DUSK) .OR.
1405
2 (DAWN .GT.DUSK .AND.
1406
3 (DJOBS(J).GT.DAWN .OR. DJOBS(J).LT.DUSK) ) ) THEN
1416
IF (DTYPE(J).EQ.'Y')THEN
1417
CALL PLOT(-1,DJOBS(J),SIGNAL(J),A1)
1420
IF (TWILIT) CARD(LWORD(CARD)+5:)='t = twilight sky'
1423
CALL PLOT(1,0.,0.,'+')
1424
WRITE(CARD,'(A,F7.0)')' TIME (MJD =',DTZERO(NIGHT)
1425
CARD(LWORD(CARD):)=' + decimal of day) -->'
1426
CALL RTNCON(CARD,LWORD(CARD))
1428
C End of TIME plot; begin AIRMASS plot.
1431
CARD(:32)='SKY brightness vs. AIRMASS, for'
1432
CARD(33:LEN+33)=BANDS(NB)
1433
CARD(LEN+34:)=DATSTR
1434
CARD(LWORD(CARD)+5:)='M = Moon up'
1439
CALL PLOT(0,XLIMS,YLIMS,'L')
1444
IF (KBND(J).NE.NB) GO TO 30
1445
IF (DTYPE(J).EQ.'Y')THEN
1447
IF ( (DAWN.LT.DUSK .AND.
1448
1 DJOBS(J).GT.DAWN .AND. DJOBS(J).LT.DUSK) .OR.
1449
2 (DAWN .GT.DUSK .AND.
1450
3 (DJOBS(J).GT.DAWN .OR. DJOBS(J).LT.DUSK) ) ) THEN
1455
ELSE IF ( (RISE.LT.SET .AND.
1456
1 DJOBS(J).GT.RISE .AND. DJOBS(J).LT.SET) .OR.
1457
2 (RISE .GT.SET .AND.
1458
3 (DJOBS(J).GT.RISE .OR. DJOBS(J).LT.SET) ) ) THEN
1459
C Moon is up. Flag in KX.
1462
C Flag for model solution.
1470
CALL PLOT(-1,AIRM(J),SIGNAL(J),A1)
1471
C Set KX(2,*) to NDIA, or 0 if no diaphragm:
1473
IF (DIAPH(J).EQ.DIANAM(NDIA)) THEN
1475
NDIA2NP(NDIA,1)=NDIA2NP(NDIA,1)+1
1479
C Here only if no dias.
1482
C IF (KDIAPHRAGM.NE.-1) CALL TV('Dia.error!')
1484
C Set KX(3,*) to J for use in later plots.
1486
C Set X(1,*) to airmass.
1488
C Set X(4,*) to nearest star signal.
1489
CALL HALO(NB,DTYPE,KDIAPHRAGM,J,JBGN,JEND, X(4,K))
1491
W(K)=1./ABS(SIGNAL(J))
1494
CALL PLOT(1,50.,-50.,' ')
1495
CALL RTNCON(' AIRMASS -->',16)
1497
C End of AIRMASS plot; begin LUNAR ELONGATION plot.
1499
IF(MOONUP.LT.2) GO TO 100
1503
CARD(:41)='SKY brightness vs. LUNAR ELONGATION, for'
1504
CARD(42:LEN+42)=BANDS(NB)
1505
CARD(LEN+43:)=DATSTR
1506
CARD(LWORD(CARD)+5:)='M = Moon up'
1511
CALL PLOT(0,XLIMS,YLIMS,'L')
1514
IF (KBND(J).NE.NB) GO TO 40
1515
IF (DTYPE(J).EQ.'Y')THEN
1517
T=(DTZERO(NIGHT)-51544.5D0+DJOBS(J))/36525.
1518
UTRAD=MOD(DJOBS(J),1.)*TWOPI
1519
TLST=STUTZ+(UTRAD*1.00273791)+ELONG
1520
C MOON uses T for lunar position, TLST for topocentric ephemeris.
1522
C MOON returns lunar coords. in RASUN, DESUN in /CSUN/.
1527
XM= -COSDEL*SIN(HAMOON)
1528
YM=COSPHI*SINDEL-SINPHI*COSDEL*COSHA
1529
ZM=SINPHI*SINDEL+COSPHI*COSDEL*COSHA
1530
IF (ZM.GT.ALTS(1)) THEN
1531
C Set X(2,*) to LUNAR airmass. (ZM is cosZ(Moon).)
1532
X(2,K)=(.0096467+ZM*(.148386+ZM*1.0024342))/
1533
1 (.000303978+ZM*(.0102963+ZM*(.149864+ZM)))
1535
C Moon below horizon.
1537
C WRITE(CARD,*)'ZM=',ZM
1538
C IF(KX(1,K).EQ.1)CALL TV(CARD)
1539
IF(KX(1,K).EQ.1) SYMB(K)='m'
1542
DELAZ=AZMOON-AZIM(J)
1543
C approximate cosz(sky) as 1/airmass.
1545
1 SQRT((1.-ZM*ZM)*(1.-1./AIRM(J)**2))*COS(DELAZ)
1546
C Set X(3,*) to lunar elongation.
1548
IF (KX(1,K).EQ.1) THEN
1549
ELONGMAX=MAX(ELONGMAX,X(3,K))
1550
ELONGMIN=MIN(ELONGMIN,X(3,K))
1552
XS(K)=ACOS(COSEL)/DEGRAD
1556
CALL PLOT(0,0.,0.,'D')
1557
CALL PLOT(K,XS,YS,SYMB)
1558
CALL PLOT(0,0.,0.,'S')
1559
CALL RTNCON(' LUNAR ELONGATION (deg.) -->',32)
1564
CALL TVN('Do you want to interpolate in TIME alone, or')
1565
CALL TVN('do you want to try fitting sky with a MODEL?')
1566
CALL TV('(Enter T for Time; M to Model the sky; or)')
1567
CALL ASKN('( R to re-display the plots.)',A1)
1570
102 IF (A1.EQ.'T') THEN
1574
C restore sky data for current night and passband.
1575
IF (KBND(J).NE.NB) GO TO 105
1576
IF (DTYPE(J).EQ.'Y')THEN
1585
110 CALL TV('Enter C to replace sky data with a Constant,')
1586
CALL TVN(' L to smooth with a Linear fit,')
1587
CALL TVN(' N to use the sky Nearest the star,')
1588
CALL TVN(' P to use the sky Preceding the star,')
1589
CALL TVN(' F to use the sky Following the star,')
1590
CALL TVN(' R to Re-display plots,')
1591
CALL TVN(' or Q to Quit.')
1594
C be lazy: fit line, but force it horizontal.
1595
CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
1600
ELSE IF (A1.EQ.'L') THEN
1602
CALL TV('Not enough data to fit a line safely.')
1603
CALL TVN('Try a constant instead.')
1606
CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
1610
C for neighboring-sample subtraction:
1613
C use NEAREST neighbor.
1614
CALL TV('CAUTION: This method is not robust!')
1615
CALL SKYNBR(NB,DTYPE,KDIAPHRAGM,JBGN,JEND)
1616
ELSE IF (A1.EQ.'P') THEN
1618
CALL TV('CAUTION: This method is not robust!')
1619
CALL SKYADJ(NB,DTYPE,KDIAPHRAGM,JBGN,JEND,NERROR)
1620
IF (NERROR.GT.0) GO TO 2
1621
ELSE IF (A1.EQ.'F') THEN
1622
C STAR comes first, so scan in reverse order.
1623
CALL TV('CAUTION: This method is not robust!')
1624
CALL SKYADJ(NB,DTYPE,KDIAPHRAGM,JEND,JBGN,NERROR)
1625
IF (NERROR.GT.0) GO TO 2
1626
ELSE IF (A1.EQ.'R') THEN
1629
CALL TV('Usually N is your best choice; but')
1630
CALL TVN('if you have been careful to take sky')
1631
CALL TVN('ALWAYS before or after every star, you')
1632
CALL TVN('can choose P for sky-then-star, or')
1633
CALL TVN(' F for star-then-sky.')
1636
C skip residual plots if sky is subtracted exactly.
1637
CALL TV('(Sky samples were subtracted exactly; so')
1638
CALL TVN(' there are no residuals to plot.)')
1642
C here when linear model was used.
1646
C PG(2) was set above.
1653
C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
1654
CALL DLSQ(YPLIN,K,2, 1, 2, NIX, 1, 0, 1.E-3)
1662
IF (KBND(J).NE.NB) GO TO 130
1663
C subtract fitted sky:
1664
SIGNAL(J)=SIGNAL(J)-((DJOBS(J)-XBAR)*SLOPE+YBAR)
1665
IF (DTYPE(J).EQ.'Y') THEN
1666
YMAX=MAX(YMAX,SIGNAL(J))
1667
YMIN=MIN(YMIN,SIGNAL(J))
1674
CARD(:26)='Sky resids. vs. TIME for'
1677
C set plot limits....
1678
XLIMS(1)=DJOBS(JMIN)
1679
XLIMS(2)=DJOBS(JMAX)
1682
CALL PLOT(0,XLIMS,YLIMS,'L')
1684
IF (KBND(J).NE.NB) GO TO 140
1685
CALL PLOT(-1,DJOBS(J),SIGNAL(J),'*')
1688
CALL PLOT(1,0.,0.,'+')
1689
CALL RTNCON(' Time (decimal days) -->',28)
1692
ELSE IF (A1.EQ.'M') THEN
1695
ELSE IF (A1.EQ.'R') THEN
1698
ELSE IF (A1.EQ.'H') THEN
1700
CALL TV('Reply M to detect bad sky measurements.')
1701
CALL TV('Reply T to use a neighboring sky sample,')
1702
CALL TVN('either before, after, or nearest the star.')
1709
ELSE IF (K.EQ.1) THEN
1710
CALL TV('Only 1 datum available for this night!')
1715
1 SIGNAL(J)=SIGNAL(J)-SIGMIN
1720
CALL TV('NO SKY DATA AVAILABLE for this night')
1721
Call ASK('Do you want to continue?',A1)
1730
C Now ***** MODEL ***** and interpolate sky.
1734
C remove const. term:
1738
C set up diaphragm factors:
1742
DO 157 K=12,11+NDIAS
1745
WRITE(CARD,*)NDIA2NP(NDIA,1),' data in dia.',DIANAM(NDIA)
1747
C Set NPARMS to largest one actually used.
1748
IF (NDIA2NP(NDIA,1).GT.0) NPARMS=K
1749
C Find leading diaphragm...
1750
IF (NDIA2NP(NDIA,1).GT.NN) THEN
1751
C flag NDIA as possible NREFDIA.
1761
IF (X(4,J).GT.0.) THEN
1767
CALL ROBLIN(XS,YS,K, XBAR, YBAR, SLOPE)
1768
C Fraction of star included:
1769
IF (SLOPE.GT.0. .AND. SLOPE.LT.0.01) THEN
1772
C ZM is max. value of XS.
1774
CALL CENTER('Sky vs. corresponding Star signal:')
1776
CALL PLOT(-K,XS,YS,' ')
1777
CALL PLOT(0,0.,0.,'O')
1780
YM=YBAR+SLOPE*(XM-XBAR)
1781
CALL PLOT(-1,XM,YM,'+')
1783
CALL PLOT(0,0.,1.,'O')
1785
IF (X(4,J).EQ.0.) GO TO 166
1786
CALL PLOT(-1,X(4,J),Y(J),SYMB(J))
1788
CALL PLOT(1,0.,0.,'+')
1789
CARD=' Star signal --> (Sky is 0.xxxx of Star)'
1790
WRITE (CARD(32:38),'(F7.4)') SLOPE
1791
CALL RTNCON(CARD,49)
1793
CALL TV('Contribution of star to sky measures not detected.')
1797
C Set default order:
1806
C Start for dark sky:
1807
C Fit rough line to dark sky:
1810
IF (KX(1,J).EQ.0) THEN
1814
YS(K)=Y(J) - PG(11)*X(4,J)
1818
C zero weights if Moon is up, or twilight:
1824
IF (NPTS-MOONUP-NTWI.GT.6 .AND. AIRMAX.GT.1.5) THEN
1827
C CALL PLOT(-NN,XS,YS,'*')
1828
C NN data now set for moonless part.
1832
C Set order of relaxation:
1837
C Fit line to data vs. airmass:
1838
CALL ROBLIN(XS,YS,NN, XBAR,YBAR,SLOPE)
1841
C CALL PLOT(-1, XS(K), YBAR+SLOPE*(XS(K)-XBAR),'-')
1843
C CALL PLOT(1,0.,0.,'+')
1844
C CALL RTNCON(' AIRMASS -->',13)
1847
IF (YBAR.LE.XBAR*SLOPE) THEN
1848
C Negative or zero intercept; form can't be fit.
1849
C CALL TV('Negative intercept.')
1854
C Positive intercept; OK.
1855
IF (SLOPE.GT.0.5*YBAR/XBAR) THEN
1857
C CALL TV('P(4) small.')
1858
PG(4)=(SQRT(YBAR/(XBAR*SLOPE))-1.)/XBAR
1862
C CALL TV('P(4) big; using P(5).')
1864
PG(5)=(YBAR/XBAR-SLOPE)/((YBAR+SLOPE*XBAR)*XBAR)
1865
PG(2)=(YBAR/XBAR+YBAR*PG(4)+YBAR*XBAR*PG(5))
1870
C Override trial line with a priori values:
1872
PG(5)=0.5*PG(4)*PG(4)
1873
PG(2)=YBAR/(1.+(PG(4)+PG(5)*XBAR)*XBAR)
1875
C Not enough data or range to fit.
1876
CALL TV('Can''t fit dark sky.')
1879
PG(5)=0.5*PG(4)*PG(4)
1881
PG(3)=PG(5)*PG(2)*0.1
1882
ZENITH=PG(1)+(PG(2)+PG(3))/(1.+PG(4)+PG(5))
1883
IF (NPTS-MOONUP-NTWI.LT.6) GO TO 200
1887
C hold all fixed but P(2):
1889
CALL FIXP(J,PG(J),NIX)
1893
C Relax successive parameters:
1895
CARD='Solving for x sky parameters...'
1896
C IF(L.EQ.3) CARD(24:24)=' '
1897
WRITE(CARD(12:13),'(I2)') L-1
1899
IF (NPTS+NIX.LT.NPARMS+4) GO TO 191
1901
C Fix all sky parms.:
1906
IF (NDIAS.GT.0) THEN
1907
C Fix reference diaphragm:
1909
CALL FIXP(K,PG(K),NIX)
1911
C WRITE(CARD,*)'Holding PG(',K,') at',PG(K)
1915
C See if we can release P(11), the star fraction:
1916
IF (NPTS.GT.2*MOONUP) CALL UNFIXP(11,NIX)
1918
C First, vary ONLY the new parameter.
1919
CALL UNFIXP(LORDER(L-1),NIX)
1921
C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
1922
CALL DLSQ(YPSKY,NPTS,NPARMS, 1, NPARMS,NIX, 4, 2, 1.E-3)
1924
C Next, see if result is acceptable:
1925
IF (WVAR.GT.0.D0) PG(LORDER(L-1))=P(LORDER(L-1))
1926
C Check each case individually:
1927
IF (P(2).LT.0.D0) CALL FIXP(2,0.,NIX)
1928
IF (P(3).LT.0.D0) CALL FIXP(3,0.,NIX)
1929
IF (P(5).LT.0.D0) CALL FIXP(5,0.,NIX)
1930
IF (P(11).LT.0.D0) CALL FIXP(11,0.,NIX)
1933
CALL UNFIXP(LORDER(K-1),NIX)
1936
C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
1937
CALL DLSQ(YPSKY,NPTS,NPARMS, 1, NPARMS,NIX, 4, 2, 1.E-3)
1939
C Give up if solution failed:
1940
IF (WVAR.EQ.0.D0 .OR. WVAR.GT.1.2*OLDVAR) GO TO 189
1944
C CALL PLOT(-NN,XS,YS,'*')
1946
C remember that YC is DP...
1947
C IF (W(K).GT.0.) CALL PLOT(-1,X(1,K),REAL(YC(K)),'C')
1949
C CALL PLOT(1,0.,0.,'+')
1950
C CALL RTNCON(' AIRMASS -->',13)
1952
C Check solution for plausibility:
1953
IF (P(2).LE.0.D0 .OR. P(3).LT.0.D0 .OR. P(5).LT.0.D0 .OR.
1954
1 (P(4).LT.0.D0 .AND. P(4)*P(4).GE.4.D0*P(5)) .OR.
1955
2 (P(4).GT.100.D0 .OR. P(5).GE.100.D0) .OR.
1956
3 (P(11).GT.0.01D0 .OR. P(11).LT.0.D0) ) THEN
1963
C Adjust weights of moonless data:
1964
ZENITH=PG(1)+(PG(2)+PG(3))/(1.+PG(4)+PG(5))
1966
IF (KX(1,K).EQ.0) W(K)=EXPTIM(KX(3,K))/(ABS(YC(K)) + ZENITH)
1968
C Update parameters:
1972
C WRITE(CARD,'(''P('',I2,'') ='',G9.2)') K,P(K)
1974
C WRITE(CARD,'(4X,''+/-'',G9.2)') SP(K)
1979
C IF (SP(11).GT.0.) THEN
1980
C WRITE(CARD,'(''P(11) ='',G9.2)') P(11)
1982
C WRITE(CARD,'(4X,''+/-'',G9.2)') SP(11)
1987
IF (P(K).GT.0.) PG(K)=P(K)
1990
WRITE(CARD,'(I5,A)') NFIT,'-parameter fit accepted.'
1993
C Release next parameter and try again.
1994
189 CALL UNFIXP(LORDER(L),NIX)
2002
WRITE (CARD,'(A,I3,A,I4,2A)')
2003
1 'Fit',NFIT,' parameters to',NN,' dark-sky data in ',BANDS(NB)
2006
ZENITH=PG(1)+(PG(2)+PG(3))/(1.+PG(4)+PG(5))
2012
YLIMS(2)=1.4*MAX(ZENITH, PG(1)+((PG(2)+AIRMAX*PG(3))*AIRMAX)/
2013
1 (1.+AIRMAX*(PG(4)+AIRMAX*PG(5))),
2014
2 PG(1)+((PG(2)+XM*PG(3))*XM)/(1.+XM*(PG(4)+XM*PG(5))))
2015
CALL PLOT(0,XLIMS,YLIMS,'L')
2016
CALL PLOT(0,0.,0.,'O')
2019
XS(1)=1.+(AIRMAX-1.)*.025*K
2020
YS(1)=PG(1)+((PG(2)+XS(1)*PG(3))*XS(1))/
2021
1 (1.+XS(1)*(PG(4)+XS(1)*PG(5)))
2022
CALL PLOT(-1,XS,YS,'+')
2024
CALL PLOT(0,0.,1.,'O')
2027
IF(KX(1,K).EQ.0) CALL PLOT(-1,X(1,K),Y(K)-PG(11)*X(4,K), '*')
2029
CALL PLOT(1,0.,0.,'+')
2030
CALL RTNCON(' Airmass -->',14)
2032
C Start for MOONLIGHT:
2034
C Rough estimates for large-airmass parameters:
2035
200 PG(9)=0.35*EXTIN(NB)
2036
PG(10)=EXP(-EXTIN(NB))
2038
C Fit rough fcn. to lunar aureole:
2041
IF (KX(1,J).EQ.1) THEN
2047
C remove dark-sky component...
2048
YS(K)=(Y(J) - (PG(1)+((PG(2)+X(1,J)*PG(3))*X(1,J))/
2049
1 (1.+X(1,J)*(PG(4)+X(1,J)*PG(5))))
2050
C ... remove star halo ...
2052
C ...and scale to unit airmass.
2053
3 ) * X(3,J)/(X(1,J)*(EXP(-PG(9)*XM)+PG(10)/ZM))
2054
C partial allowance for photon noise:
2055
W(J)=EXPTIM(KX(3,J))/(ABS(Y(J))+ZENITH)
2057
C Keep old wt. for dark sky.
2061
IF (MOONUP.GT.6 .AND. ELONGMAX-ELONGMIN.GT.0.2) THEN
2062
C Sort on elongation & use bottom of list for aureole:
2065
C Select points in aureole:
2066
IF (XS(J).GT.1.) GO TO 212
2068
C Whole set is in aureole. Use 1st half:
2070
212 IF (J.EQ.1) THEN
2071
C No aureole data. Pick something:
2075
C WRITE(CARD,*) J,' points in initial fit.'
2077
C First J now suitable to estimate PG(6).
2078
CALL ROBLIN(XS,YS,J, XBAR,YBAR,SLOPE)
2079
PG(6)=YBAR-XBAR*SLOPE
2081
IF (PG(6).LT.0.D0) THEN
2082
C Aureole not detectable.
2086
C Now move down *last* third...
2089
C ...and divide by elongation.
2090
YS(J)=YS(K-K/3+J)/XS(J)
2092
CALL ROBLIN(XS,YS,K/3, XBAR,YBAR,SLOPE)
2095
C Not enough data to fit with Moon.
2099
PG(7)=(YS((K+1)/2)+YS(K/2+1))*0.5
2107
C PG's now set for moonlight terms.
2109
IF (ELONGMAX-ELONGMIN.LT.2.) THEN
2110
C Rearrange order of solving:
2111
IF (ELONGMIN.GT.1.5) THEN
2115
IF (ELONGMIN.GT.2.) THEN
2119
ELSE IF (ELONGMIN.GT.0.5) THEN
2122
IF (ELONGMAX.LT.1.) THEN
2128
ELSE IF (ELONGMAX.LT.1.0) THEN
2138
IF (MOONUP.LT.3) GO TO 241
2140
C Relax parameters one by one:
2145
IF (NDIAS.GT.0) THEN
2146
C Fix reference diaphragm:
2148
CALL FIXP(K,PG(K),NIX)
2151
C Release P(11), the star fraction:
2154
CARD='Solving for x moonlight parameters...'
2155
IF(J.EQ.6) CARD(24:24)=' '
2156
WRITE(CARD(12:13),'(I2)') J-5
2159
C First, vary ONLY the new parameter.
2160
CALL UNFIXP(LORDER(J),NIX)
2161
IF (MOONUP.LT.NPARMS+4-NIX) GO TO 241
2163
C subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
2164
CALL DLSQ(YPSKY,NPTS,NPARMS, 1,NPARMS,NIX, 4, 2, 1.E-3)
2166
C Next, see if result is acceptable:
2167
IF (WVAR.GT.0.D0) PG(LORDER(J))=P(LORDER(J))
2168
C Check each case individually:
2169
IF (P(6).LT.0.D0) CALL FIXP(6,0.,NIX)
2170
IF (P(8).LT.0.D0) CALL FIXP(8,0.,NIX)
2171
IF (P(9).LT.0.D0) CALL FIXP(9,0.,NIX)
2172
IF (P(10).LT.0.D0) CALL FIXP(10,0.,NIX)
2173
IF (P(9).GT.9.D0 .AND. P(10).EQ.0.D0) CALL FIXP(9,9.,NIX)
2174
IF (P(9).GT.10.D0) CALL FIXP(9,1.,NIX)
2175
IF (P(11).LT.0.D0) CALL FIXP(11,0.,NIX)
2177
C Then, vary parameters 6 to J.
2179
CALL UNFIXP(LORDER(K),NIX)
2181
IF (MOONUP.LT.NPARMS+4-NIX) GO TO 241
2183
C subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
2184
CALL DLSQ(YPSKY,NPTS,NPARMS, 1,NPARMS,NIX, 4, 2, 1.E-3)
2186
C Re-extract data for plots:
2191
IF (KX(1,L).EQ.1) THEN
2197
C remove dark-sky component...
2198
YS(K)=(Y(L) - (P(1)+((P(2)+X(1,L)*P(3))*X(1,L))/
2199
1 (1.+X(1,L)*(P(4)+X(1,L)*P(5))))
2200
C ...subtract star halo, and scale to unit airmass.
2201
2 - PG(11)*X(4,L) ) / (X(1,L)*(EXP(-P(9)*XM) + P(10)/ZM))
2202
IF(NDIAS.GT.1) YS(K)=YS(K)/P(KX(2,L)+11)
2203
YMAX=MAX(YMAX,YS(K))
2204
YMIN=MIN(YMIN,YS(K))
2211
IF (NIX.GT.0) GO TO 231
2213
C Fix up plot limits:
2217
C Adjust for possible outliers:
2219
C (Note reversal of Y and X.)
2220
YM=(YS((K+1)/2)+YS(K/2+1))*0.5
2224
C ZM is interquartile range.
2225
YMAX=MIN(YMAX,YM+4.*ZM)
2226
YMIN=MAX(YMIN,YM-4.*ZM)
2233
WRITE(CARD,'(I5,A,I4,1X,2A)')
2234
1 NPARMS-NIX,'-parameter fit ( + + + ) to',MOONUP,
2235
2 BANDS(NB)(:LWORD(BANDS(NB))),' data ( * ) with Moon'
2237
CALL TVN('Normalized Sky')
2238
CALL PLOT(0,0.,0.,'O')
2240
CALL PLOT(0,XLIMS,YLIMS,'L')
2241
IF (YLIMS(1).LT.0.) CALL XAXIS(XLIMS)
2242
C Plot fit to aureole terms:
2244
XM=0.1+(.02+.003*K)*K
2245
YM=(P(6)/XM + P(7) + P(8)*XM)
2246
CALL PLOT(-1,XM,YM,'+')
2248
CALL PLOT(0,0.,1.,'O')
2249
CALL PLOT(NN,XS,YS,'*')
2250
CALL RTNCON(' LUNAR ELONGATION (rad.) -->',32)
2252
C Give up if solution failed:
2253
231 IF (WVAR.EQ.0.D0 .OR. WVAR.GT.1.2*OLDVAR) GO TO 240
2254
C Tests for plausibility:
2255
IF (P(6).LT.0.D0 .OR. P(8).LT.0.D0 .OR.
2256
1 P(9).LT.0.D0 .OR. P(10).LT.0.D0 .OR.
2257
2 (P(9).GT.10.D0 .AND. P(10).EQ.0.D0) .OR.
2258
3 (P(11).GT.0.01D0 .OR. P(11).LT.0.D0) ) GO TO 240
2264
C partial allowance for photon noise:
2265
IF(KX(1,K).EQ.1) W(K)=EXPTIM(KX(3,K))/(ABS(Y(K))+ZENITH)
2267
C Update parameters:
2271
C WRITE(CARD,'(''P('',I2,'') ='',G9.2)') K,P(K)
2273
C WRITE(CARD,'(4X,''+/-'',G9.2)') SP(K)
2276
C Wings of star image:
2277
IF (P(11).GT.0.D0) PG(11)=P(11)
2278
IF (SP(11).GT.0.) THEN
2279
C WRITE(CARD,'(''P(11) ='',G9.2)') P(11)
2281
C WRITE(CARD,'(4X,''+/-'',G9.2)') SP(11)
2286
IF (P(K).GT.0.) PG(K)=P(K)
2290
WRITE(CARD,'(I5,A)') NFIT,'-parameter fit accepted.'
2296
C Re-extract data for final plot:
2301
IF (KX(1,L).EQ.1) THEN
2307
C remove dark-sky component...
2308
YS(K)=(Y(L) - (PG(1)+((PG(2)+X(1,L)*PG(3))*X(1,L))/
2309
1 (1.+X(1,L)*(PG(4)+X(1,L)*PG(5))))
2310
C ...subtract star halo, and scale to unit airmass.
2311
2 - PG(11)*X(4,L) ) / (X(1,L)*(EXP(-PG(9)*XM) + PG(10)/ZM))
2312
YMAX=MAX(YMAX,YS(K))
2313
YMIN=MIN(YMIN,YS(K))
2317
C Fix up plot limits:
2320
IF (K.EQ.0) GO TO 300
2321
C Adjust for possible outliers:
2323
C (Note reversal of Y and X.)
2324
YM=(YS((K+1)/2)+YS(K/2+1))*0.5
2328
C ZM is interquartile range.
2329
YMAX=MIN(YMAX,YM+3.*ZM)
2330
YMIN=MAX(YMIN,YM-3.*ZM)
2331
IF (YMAX.EQ.YMIN) GO TO 300
2338
WRITE(CARD,'(I5,A,I4,1X,2A)')
2339
1 NFIT,'-parameter fit ( + + + ) to',MOONUP,
2340
2 BANDS(NB)(:LWORD(BANDS(NB))),' data ( * ) with Moon'
2342
CALL TVN('Normalized Sky')
2343
CALL PLOT(0,0.,0.,'O')
2345
CALL PLOT(0,XLIMS,YLIMS,'L')
2346
IF (YLIMS(1).LT.0.) CALL XAXIS(XLIMS)
2347
C Plot fit to aureole terms:
2349
XM=ELONGMIN+(ELONGMAX-ELONGMIN)*(.02+.003*K)*K
2350
YM=(PG(6)/XM + PG(7) + PG(8)*XM)
2351
CALL PLOT(-1,XM,YM,'+')
2353
CALL PLOT(0,0.,1.,'O')
2354
CALL PLOT(NN,XS,YS,'*')
2355
CALL RTNCON(' LUNAR ELONGATION (rad.) -->',32)
2360
C Now display fitted equation:
2365
CALL SKYEQN(PG,SP,CARD)
2368
C Don't forget diaphragm scaling:
2370
IF (NPARMS.GT.11) THEN
2372
CALL TV('Diaphragm scaling:')
2375
WRITE(CARD,'(A,F7.3,'' +/-'',F7.3)') DIANAM(J-11),P(J),SP(J)
2380
C Now display overall fit:
2382
C Hold all fixed; evaluate YC's:
2388
C subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
2389
CALL DLSQ(YPSKY,NPTS,NPARMS, 1, NPARMS,NIX, 4, 2, 1.E-1)
2393
CARD='Measured vs. model ( + + ) sky brightness in '//BANDS(NB)
2394
CARD(LWORD(CARD)+1:)=':'
2396
CALL TVN('Observed sky')
2403
IF (KX(1,J).LT.0) GO TO 360
2407
XLIMS(1)=MIN(XLIMS(1),XS(K))
2408
IF (KX(1,J).EQ.0) THEN
2410
ELSE IF (KX(1,J).EQ.1) THEN
2412
C ELSE IF (KX(1,J).EQ.-1) THEN
2419
IF (XS(K).LT.2.*XS(K-1)) THEN
2424
YLIMS(1)=0.8*XLIMS(1)
2425
YLIMS(2)=1.2*XLIMS(2)
2426
CALL PLOT(0,XLIMS,YLIMS,'L')
2427
C Draw unit diagonal:
2430
CALL PLOT(-1,XM,XM,'+')
2432
CALL PLOT(0,0.,0.,'D')
2433
CALL PLOT(K,XS,YS,SYMB)
2434
CALL RTNCON(' Calculated sky -->',23)
2436
C Plot star halo resids.:
2440
CARD='Residuals vs. Star signal in '//BANDS(NB)
2441
CARD(LWORD(CARD)+1:)=':'
2443
CALL TVN(' Residual')
2450
IF (KX(1,J).LT.0) GO TO 380
2454
XLIMS(2)=MAX(XLIMS(2),XS(K))
2455
YLIMS(1)=MIN(YLIMS(1),YS(K))
2456
YLIMS(2)=MAX(YLIMS(2),YS(K))
2458
CALL PLOT(0,XLIMS,YLIMS,'L')
2460
CALL PLOT(K,XS,YS,SYMB)
2461
CALL PLOT(0,0.,0.,'S')
2462
CARD=' Star signal --> (Sky is 0.xxxx of Star)'
2463
WRITE (CARD(32:38),'(F7.4)') PG(11)
2464
CALL RTNCON(CARD,49)
2466
C Plot resids. vs. time:
2470
CARD='O/C Ratio vs. Time in '//BANDS(NB)
2471
CARD(LWORD(CARD)+1:)=':'
2474
XLIMS(1)=DJOBS(JDUSK)
2475
XLIMS(2)=DJOBS(JDAWN)
2478
CALL PLOT(0,XLIMS,YLIMS,'L')
2479
C Draw "axis" at 1.0:
2481
XM=XLIMS(1)+(XLIMS(2)-XLIMS(1))*J/80.
2482
CALL PLOT(-1,XM,1.,'-')
2485
IF (KX(1,J).LT.0) GO TO 390
2486
IF (YC(J).EQ.0.D0) YC(J)=1.E-6
2487
CALL PLOT(-1,DJOBS(KX(3,J)),Y(J)/REAL(YC(J)),SYMB(J))
2490
CALL PLOT(1,10.,0.,'+')
2491
WRITE(CARD,'(A,F7.0)')' TIME (MJD =',DTZERO(NIGHT)
2492
CARD(LWORD(CARD):)=' + decimal of day) -->'
2493
CALL RTNCON(CARD,LWORD(CARD))
2500
CARD='Sky anomalies in '//BANDS(NB)
2501
WRITE(CARD(LWORD(CARD)+2:),'(''on MJD'',F7.0)')
2502
1 DTZERO(NIGHT)+INT(DJOBS((JBGN+JEND)/2))
2504
CARD='DECIMAL SKY LUNAR LUNAR STELLAR'
2505
CARD(52:)='SKY SKY O/C'
2507
CARD=' DAY AIRMASS AIRMASS ELONGATION SIGNAL OBSERVED'
2508
CARD(59:)='CALCULATED RATIO'
2510
WRITE(CARD,'(''------'',7(3X,''-------''))')
2514
IF (KX(1,J).GE.0 .AND.
2515
1 W(J)*(Y(J)-REAL(YC(J)))**2 .GT. SCALE*36.) THEN
2516
WRITE(CARD,'(F6.4,F9.3,2F10.3,1PE12.2,0P2F10.3,F8.3)')
2517
1 DJOBS(KX(3,J)),(X(K,J),K=1,4),Y(J),YC(J),Y(J)/YC(J)
2518
IF (KX(1,J).NE.1) THEN
2528
WRITE(CARD,'(I3,'' anomalous sky data'')') NN
2529
IF (NN.EQ.1) CARD(22:23)='um'
2537
1 'Do you want to subtract the fitted model from stellar data?'
2539
IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN
2541
ELSE IF (A1.EQ.'N') THEN
2542
C Take short-cut to 110.
2545
ELSE IF (A1.EQ.'R') THEN
2547
ELSE IF (A1.EQ.'H') THEN
2548
IF (NPTS.GT.20*NN .AND. SCALE.LT.0.1*ZENITH) THEN
2549
C Clearly a good fit.
2550
CALL TV(' YES, that''s a good fit.')
2551
ELSE IF (NPTS.LT.10*NN .OR. SCALE.GT.0.2*ZENITH) THEN
2552
C Clearly a bad fit.
2553
CALL TV(' NO, that''s a POOR fit.')
2557
WRITE(CARD,'(I5,''% rejects; typical residual is'',F5.2)')
2558
1 (NN*100+NPTS/2)/NPTS,SCALE/ZENITH
2559
CARD(LWORD(CARD)+2:)='of zenith brightness.'
2561
CALL TV('Hard to say. You may want to try it both ways,')
2562
CALL TVN('and compare the results.')
2567
CALL TV('Enter either YES or NO, or')
2568
CALL TVN('H for Help, or R to Re-display plots.')
2573
C Now subtract sky fit from data:
2575
C JBGN, JEND index the first and last data of ANY kind for the night.
2576
C JMIN, JMAX index the first and last SKY data for the night & band.
2580
CALL TV('Using nearest-neighbor subtraction in twilight...')
2581
CALL SKYNBR(NB,DTYPE,KDIAPHRAGM,JBGN,JDUSK)
2582
CALL SKYNBR(NB,DTYPE,KDIAPHRAGM,JDAWN,JEND)
2585
C Subtract fitted fcn. from rest:
2586
CALL TV('Using fitted model...')
2587
DO 600 J=JDUSK,JDAWN
2588
IF (NB.NE.KBND(J)) GO TO 600
2590
IF ( (RISE.LT.SET .AND.
2591
1 DJOBS(J).GT.RISE .AND. DJOBS(J).LT.SET) .OR.
2592
2 (RISE .GT.SET .AND.
2593
3 (DJOBS(J).GT.RISE .OR. DJOBS(J).LT.SET) ) ) THEN
2596
T=(DTZERO(NIGHT)-51544.5D0+DJOBS(J))/36525.
2597
UTRAD=MOD(DJOBS(J),1.)*TWOPI
2598
TLST=STUTZ+(UTRAD*1.00273791)+ELONG
2599
C MOON uses T for lunar position, TLST for topocentric ephemeris.
2601
C MOON returns lunar coords. in RASUN, DESUN in /CSUN/.
2606
XM= -COSDEL*SIN(HAMOON)
2607
YM=COSPHI*SINDEL-SINPHI*COSDEL*COSHA
2608
ZM=SINPHI*SINDEL+COSPHI*COSDEL*COSHA
2610
C Set Z(2) to LUNAR airmass.
2613
C Moon below horizon.
2617
DELAZ=AZMOON-AZIM(J)
2618
C approximate cosz(sky) as 1/airmass.
2620
1 SQRT((1.-ZM*ZM)*(1.-1./AIRM(J)**2))*COS(DELAZ)
2621
C Set Z(3) to lunar elongation.
2629
IF (DTYPE(J).EQ.'S') THEN
2632
ELSE IF (DTYPE(J).EQ.'Y') THEN
2634
CALL HALO(NB,DTYPE,KDIAPHRAGM,J,JBGN,JEND, Z(4))
2640
C Put NDIA in KZ(2):
2642
IF (DIAPH(J).EQ.DIANAM(NDIA)) THEN
2650
C Z's and KZ(1) are now set for YPSKY:
2652
C Sky brightness now returned in YT.
2653
SIGNAL(J)=SIGNAL(J)-YT
2657
C Prepare for next night:
2662
C END of loops over bands and nights.
2668
C Fits Y = [Y1 + Y2 + Y3 + P(1)] * P(diaphragm), where
2672
C P(2)*Z(1) + P(3)*Z(1)
2673
C Y1 = --------------------------- is the sky background, and
2675
C 1 + P(4)*Z(1) + P(5)*Z(1)
2679
C Y2 = Z(1) * [P(6)/Z(3) + P(7) + P(8)*Z(3)] * {exp[-P(9)*v1] + P(10)/v3}
2681
C is the scattered moonlight.
2684
C Y3 = Z(4) * P(11) is the part of the star included in sky measures.
2688
C P(1) is a possible zero-offset, if dark was not subtracted.
2691
C Independent variables:
2692
C ---------------------
2694
C Z(1) = observed airmass v1 = Z(1) + Z(2) = sum of airmasses
2696
C Z(2) = lunar airmass v3 = Z(1)*Z(2) = product of airmasses
2698
C Z(3) = lunar elongation
2700
C Z(4) = nearest star signal, or zero if none found.
2703
C KZ(1) = Moon flag (1 if Moon up)
2705
C KZ(2) = diaphragm number (KDIA)
2712
C Includes for dlsq:
2713
INCLUDE 'MID_REL_INCL:kingfit.inc'
2714
INCLUDE 'MID_REL_INCL:dlsq.inc'
2716
DOUBLE PRECISION DEN, V1, V3, FIRST, SECOND, FZ1, EXPFAC, Y2
2720
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
2723
IF (KZ(1).LT.0) THEN
2724
C Twilight: Not in model.
2732
C P(2)*Z(1) + P(3)*Z(1)
2733
C Y1 = --------------------------- is the sky background...
2735
C 1 + P(4)*Z(1) + P(5)*Z(1)
2746
DEN = 1.D0 + (P(4) + P(5)*Z(1))*Z(1)
2747
YT = (P(2) + P(3)*Z(1))*Z(1)/DEN
2751
PART(3)=Z(1)*PART(2)
2752
PART(4)= -YT*PART(2)
2753
PART(5)= -YT*PART(3)
2756
YT = YT + P(1) + P(11)*Z(4)
2758
IF (KZ(1).EQ.0) THEN
2759
C Don't need moonlight, so finish up:
2760
IF (KZ(2).GT.0) THEN
2761
C Scale old partials by diaphragm factor...
2764
PART(IP(I))=PART(IP(I))*P(KDIA)
2766
C Scale by diaphragm factor, and add to KIP:
2776
C Here if we need to add moonlight:
2779
C Y2 = Z(1) * [P(6)/Z(3) + P(7) + P(8)*Z(3)] * {exp[-P(9)*v1] + P(10)/v3}
2781
C is the scattered moonlight.
2793
FIRST=P(6)/Z(3) + P(7) + P(8)*Z(3)
2795
IF (-P(9)*V1.GT.30.D0) THEN
2796
C try to prevent floating-point exceptions.
2799
EXPFAC=EXP(-P(9)*V1)
2801
SECOND=EXPFAC + P(10)/V3
2806
PART(6)=PART(7)/Z(3)
2807
PART(8)=Z(3)*PART(7)
2808
PART(9)= -V1*EXPFAC*FZ1
2813
C Diaphragm scaling:
2814
IF (KZ(2).GT.0) THEN
2815
C Scale old partials by diaphragm factor...
2818
PART(IP(I))=PART(IP(I))*P(KDIA)
2820
C Scale by diaphragm factor, and add to KIP:
2829
SUBROUTINE SKYADJ(NB,DTYPE,KDIAPHRAGM,JBGN,JEND,NERROR)
2831
C Picks an ADJACENT sky sample, and subtracts it from star.
2836
REAL SKY, PHNOISE, ERREST
2837
INTEGER NB, KDIAPHRAGM, JBGN, JEND, NERROR, J3, J,
2840
INCLUDE 'MID_REL_INCL:obs.inc'
2842
INCLUDE 'MID_REL_INCL:mstars.inc'
2843
C PARAMETER (MSTARS=1650)
2844
C MSTARS is max.number of star-catalog entries.
2845
C commons for star catalog:
2847
COMMON /SCATA/ STARS(MSTARS)
2849
CHARACTER DTYPE(MXOBS)
2851
CHARACTER A, CARD*48, PRECEDED*8, FOLLOWING*9, DSKY*4
2858
C Decide which way to scan:
2860
IF (JBGN.LT.JEND) THEN
2863
FOLLOWING='FOLLOWING'
2867
FOLLOWING='PRECEDING'
2870
DO 134 J=JBGN,JEND,J3
2872
IF (KBND(J).NE.NB .OR. NSTR(J).LT.0) GO TO 134
2873
C here if bands match.
2874
IF (DTYPE(J).EQ.'Y') THEN
2875
C this is sky; reset it.
2879
IF (KDIAPHRAGM.NE.-1) THEN
2880
C flag sky diaphragm code; expect star to match.
2884
C this is star; subtract sky.
2886
IF (SKY.EQ.-3.E33) THEN
2887
C ERROR: no previous sky!
2890
CARD(LWORD(CARD)+2:LWORD(CARD)+6)='at'
2891
WRITE(CARD(LWORD(CARD)+2:),'(F7.4)')DJOBS(J)
2893
CARD='!! Star not '//PRECEDED//' by sky !!'
2895
CARD='Do you want to use the '//FOLLOWING//' sky?'
2898
IF (A.EQ.'Y' .OR. A.EQ.'O')THEN
2899
C look for following sky:
2900
DO 132 JJ=J+J3,JEND,J3
2901
IF (KBND(JJ).NE.NB .OR.
2902
1 NSTR(JJ).LT.0) GO TO 132
2903
IF(DTYPE(JJ).EQ.'Y') THEN
2904
IF(KDIAPHRAGM.NE.-1 .AND.
2905
1 DIAPH(JJ).NE.DIAPH(J)) GOTO 132
2906
C fix this one, but do not set SKY.
2907
SIGNAL(J)=SIGNAL(J)-SIGNAL(JJ)
2908
C accumulate photon variances.
2909
PHVAR(J)=PHVAR(J)+PHVAR(JJ)
2911
1 ESTER(J)=SQRT(ESTER(J)**2 + ESTER(JJ)**2)
2915
CALL TV('No sky found! Quitting.')
2922
C probably normal case. Check DIAPH:
2923
IF (KDIAPHRAGM.NE.-1 .AND.
2924
1 DIAPH(J).NE.DSKY) THEN
2927
CARD(LWORD(CARD)+2:LWORD(CARD)+6)='at'
2928
WRITE(CARD(LWORD(CARD)+2:),
2932
1 'Sky and Star taken through DIFFERENT diaphragms.')
2933
CALL TVN('Sky subtraction impossible.')
2936
C NORMAL case. Subtract sky:
2937
SIGNAL(J)=SIGNAL(J)-SKY
2938
PHVAR(J)=PHVAR(J)+PHNOISE
2939
IF(ESTER(J).NE.0.) ESTER(J)=SQRT(ESTER(J)**2 + ERREST**2)
2953
999 CALL ASK('Do you want to try another sky-subtraction method?',A)
2954
IF (A.EQ.'Y' .OR. A.EQ.'O') THEN
2957
ELSE IF (A.EQ.'H') THEN
2958
CALL TV('The photometric system changes with')
2959
CALL TVN('diaphragm size, so you cannot combine')
2960
CALL TVN('data taken through different holes.')
2961
CALL RTNCON('Try the "Nearest" option.',25)
2965
CALL TV('Quitting.')
2970
SUBROUTINE HALO(NB,DTYPE,KDIAPHRAGM,J,JBGN,JEND, SKYUSE)
2972
C Picks the NEAREST STAR sample to sky(J), and returns it in SKYUSE.
2977
C DTYPE is 'Y' for sky, 'S' for star
2978
C KDIAPHRAGM is column pointer for diaphragm data
2979
C JBGN is observation number for start of night
2980
C JEND is observation number for end of night
2986
INTEGER NB, KDIAPHRAGM, JBGN, JEND
2988
INCLUDE 'MID_REL_INCL:obs.inc'
2990
CHARACTER DTYPE(MXOBS)
2994
REAL SKY, SKY2, SKYAIR, SKYAZ, SKYTIM, PI,TWOPI,
2995
1 SKYAIR2, SKYAZ2, SKYTIM2, S1, S2, SKYUSE
2999
C Meaning of local variables:
3001
C Previous Following
3002
C -------- ---------
3004
C Airmass SKYAIR SKYAIR2
3005
C Azimuth SKYAZ SKYAZ2
3006
C Time SKYTIM SKYTIM2
3007
C Diaphragm DSKY DSKY2
3010
INCLUDE 'MID_REL_INCL:mstars.inc'
3011
C PARAMETER (MSTARS=1650)
3012
C MSTARS is max.number of star-catalog entries.
3013
C commons for star catalog:
3015
COMMON /SCATA/ STARS(MSTARS)
3017
CHARACTER DSKY*4, DSKY2*4
3019
DATA PI/3.1415926535/
3033
C J is current sky datum in band NB.
3035
IF (DTYPE(J).EQ.'Y') THEN
3037
IF (KDIAPHRAGM.NE.-1) THEN
3038
C save new diaphragm.
3041
C no diaphragm info.
3044
CALL TV('Called HALO with non-sky datum.')
3048
C This is sky; look for PREVIOUS star.
3049
20 IF (SKY.LT.0.) THEN
3050
C No previous star. Find one.
3051
DO 28 JJ=J-1,JBGN,-1
3052
IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 28
3053
IF(DTYPE(JJ).EQ.'S') THEN
3054
IF(KDIAPHRAGM.NE.-1 .AND.
3055
1 DIAPH(JJ).NE.DIAPH(J)) GOTO 28
3056
C OK; JJ is the previous star we seek.
3069
C probably normal case. Check DIAPH:
3070
IF (KDIAPHRAGM.NE.-1 .AND. DIAPH(J).NE.DSKY) THEN
3071
C Mismatch; force search for previous star.
3075
C NORMAL case. Existing star OK.
3079
C SKY is now valid previous value. Look for FOLLOWING value:
3081
40 IF (SKY2.LT.0.) THEN
3082
C Look for next star.
3084
IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 48
3085
IF(DTYPE(JJ).EQ.'S') THEN
3086
IF(KDIAPHRAGM.NE.-1 .AND.
3087
1 DIAPH(JJ).NE.DIAPH(J)) GOTO 48
3088
C OK; JJ is the following star we seek.
3097
C No following star.
3101
C existing sky2 might still be OK. Check it:
3102
IF (SKYTIM2.GT.DJOBS(J)) THEN
3103
C Time is OK. Check diaphragm:
3104
IF (KDIAPHRAGM.NE.-1) THEN
3105
IF (DIAPH(J).EQ.DSKY2) THEN
3113
C No problem. Go on.
3116
C We are past the old second star. Move it back & get new one:
3127
C Now compare to find "nearer" star:
3129
60 IF (SKY.GT.0. .AND. SKY2.GT.0.) THEN
3130
C Everything is hunky-dory. Find closer sample:
3131
S1=ABS(SKYTIM-DJOBS(J))*20. + ABS(SKYAIR-AIRM(J))
3132
IF (ABS(SKYAZ-AZIM(J)).GT.PI) THEN
3133
S1=S1+ (SKYAIR+AIRM(J))*(TWOPI-ABS(SKYAZ-AZIM(J)))
3135
S1=S1+ (SKYAIR+AIRM(J))*ABS(SKYAZ-AZIM(J))
3137
S2=ABS(SKYTIM2-DJOBS(J))*20. + ABS(SKYAIR2-AIRM(J))
3138
IF (ABS(SKYAZ2-AZIM(J)).GT.PI) THEN
3139
S2=S2+ (SKYAIR2+AIRM(J))*(TWOPI-ABS(SKYAZ2-AZIM(J)))
3141
S2=S2+ (SKYAIR2+AIRM(J))*ABS(SKYAZ2-AZIM(J))
3156
ELSE IF (SKY.GT.0.) THEN
3157
C SKY is the only good one. Use it.
3159
ELSE IF (SKY2.GT.0.) THEN
3160
C SKY2 is the only good one. Use it.
3163
C NEITHER star is any good.
3164
C set value to zero and return:
3169
C OK, we have now rigged SKYUSE to be the value to use. Report it.
3177
SUBROUTINE SKYNBR(NB,DTYPE,KDIAPHRAGM,JBGN,JEND)
3179
C Picks the NEAREST sky sample, and subtracts it from star.
3184
C DTYPE is 'Y' for sky
3185
C KDIAPHRAGM is column pointer for diaphragm data
3186
C JBGN is observation number for start of night
3187
C JEND is observation number for end of night
3193
INTEGER NB, KDIAPHRAGM, JBGN, JEND
3195
INCLUDE 'MID_REL_INCL:obs.inc'
3197
CHARACTER DTYPE(MXOBS)
3201
REAL SKY, SKY2, SKYAIR, SKYAZ, SKYTIM, PI,TWOPI,
3202
1 PHNOISE, ERREST, SKYAIR2, SKYAZ2, SKYTIM2,
3203
2 PHN2, ERR2, S1, S2, SKYUSE, PHNUSE, ERRUSE
3207
C Meaning of local variables:
3209
C Previous Following
3210
C -------- ---------
3212
C Airmass SKYAIR SKYAIR2
3213
C Azimuth SKYAZ SKYAZ2
3214
C Time SKYTIM SKYTIM2
3215
C Phot.noisePHNOISE PHN2
3216
C Est.error ERREST ERR2
3217
C Diaphragm DSKY DSKY2
3220
INCLUDE 'MID_REL_INCL:mstars.inc'
3221
C PARAMETER (MSTARS=1650)
3222
C MSTARS is max.number of star-catalog entries.
3223
C commons for star catalog:
3225
COMMON /SCATA/ STARS(MSTARS)
3227
CHARACTER CARD*48, DSKY*4, DSKY2*4
3229
DATA PI/3.1415926535/
3240
IF (KBND(J).NE.NB .OR. NSTR(J).LT.0) GO TO 100
3241
C here if bands match.
3243
IF (DTYPE(J).EQ.'Y') THEN
3245
IF (KDIAPHRAGM.NE.-1) THEN
3246
C save new diaphragm.
3249
C no diaphragm info.
3259
C this is star; look for PREVIOUS sky.
3260
20 IF (SKY.LT.0.) THEN
3261
C No previous sky. Find one.
3262
DO 28 JJ=J-1,JBGN,-1
3263
IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 28
3264
IF(DTYPE(JJ).EQ.'Y') THEN
3265
IF(KDIAPHRAGM.NE.-1 .AND.
3266
1 DIAPH(JJ).NE.DIAPH(J)) GOTO 28
3267
C OK; JJ is the previous sky we seek.
3282
C probably normal case. Check DIAPH:
3283
IF (KDIAPHRAGM.NE.-1 .AND. DIAPH(J).NE.DSKY) THEN
3284
C Mismatch; force search for previous sky.
3288
C NORMAL case. Existing sky OK.
3293
C SKY is now valid previous value. Look for FOLLOWING value:
3295
40 IF (SKY2.LT.0.) THEN
3296
C Look for next sky.
3298
IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 48
3299
IF(DTYPE(JJ).EQ.'Y') THEN
3300
IF(KDIAPHRAGM.NE.-1 .AND.
3301
1 DIAPH(JJ).NE.DIAPH(J)) GOTO 48
3302
C OK; JJ is the following sky we seek.
3317
C existing sky2 might still be OK. Check it:
3318
IF (SKYTIM2.GT.DJOBS(J)) THEN
3319
C Time is OK. Check diaphragm:
3320
IF (KDIAPHRAGM.NE.-1) THEN
3321
IF (DIAPH(J).EQ.DSKY2) THEN
3329
C No problem. Go on.
3332
C We are past the old second sky. Move it back & get new one:
3345
C Now compare to find "nearer" sky:
3347
60 IF (SKY.GT.0. .AND. SKY2.GT.0.) THEN
3348
C Everything is hunky-dory. Find closer sample:
3349
S1=ABS(SKYTIM-DJOBS(J))*20. + ABS(SKYAIR-AIRM(J))
3350
IF (ABS(SKYAZ-AZIM(J)).GT.PI) THEN
3351
S1=S1+ (SKYAIR+AIRM(J))*(TWOPI-ABS(SKYAZ-AZIM(J)))
3353
S1=S1+ (SKYAIR+AIRM(J))*ABS(SKYAZ-AZIM(J))
3355
S2=ABS(SKYTIM2-DJOBS(J))*20. + ABS(SKYAIR2-AIRM(J))
3356
IF (ABS(SKYAZ2-AZIM(J)).GT.PI) THEN
3357
S2=S2+ (SKYAIR2+AIRM(J))*(TWOPI-ABS(SKYAZ2-AZIM(J)))
3359
S2=S2+ (SKYAIR2+AIRM(J))*ABS(SKYAZ2-AZIM(J))
3378
ELSE IF (SKY.GT.0.) THEN
3379
C SKY is the only good one. Use it.
3383
ELSE IF (SKY2.GT.0.) THEN
3384
C SKY2 is the only good one. Use it.
3389
C AARRRGHHH!! NEITHER sky is any good!
3390
CARD='Can''t find a sky sample to match '//STARS(NSTR(J))
3392
WRITE(CARD,'(''at '',F7.4)') DJOBS(J)
3393
IF (KDIAPHRAGM.NE.-1) CARD(30:)='in diaphragm '//DIAPH(J)
3395
CALL TV('This value will be REJECTED.')
3397
C set value to zero and continue:
3403
C OK, we have now rigged SKYUSE to be the value to subtract. Do it.
3405
SIGNAL(J)=SIGNAL(J)-SKYUSE
3406
PHVAR(J)=PHVAR(J)+PHNUSE
3407
IF (ESTER(J).NE.0.) ESTER(J)=SQRT(ESTER(J)**2 +ERRUSE**2)
3416
SUBROUTINE REDSUB(BANDS,KMIN,KMAX,LASTK,MAXK,DTYPE)
3418
C Subtracts red-leak data from corresponding bands.
3419
C Picks the NEAREST red-leak sample, and subtracts it from star.
3424
REAL RLFACT, RL, RLVAR, ERVAR
3425
INTEGER KMIN, KMAX, LASTK, MAXK, NBRL, K, KRL, JBGN, NIGHT,
3426
1 KK, J, JMIN, JMAX, JEND, ISTR, JOFF, J1, J2
3428
INCLUDE 'MID_REL_INCL:obs.inc'
3430
INCLUDE 'MID_REL_INCL:mstars.inc'
3431
C PARAMETER (MSTARS=1650)
3432
C MSTARS is max.number of star-catalog entries.
3433
C commons for star catalog:
3435
COMMON /SCATA/ STARS(MSTARS)
3437
CHARACTER DTYPE(MXOBS)
3439
INCLUDE 'MID_REL_INCL:mbands.inc'
3440
C PARAMETER (MBANDS=9)
3442
PARAMETER (MXCOLR=3*MBANDS)
3443
CHARACTER *8 BANDS(MXCOLR)
3445
C Common for red leaks:
3447
COMMON /REDLK/ NBRL(MXCOLR), RLFACT(MXCOLR)
3450
C Useful bands run from KMIN to LASTK.
3451
C Measured bands run from KMIN to KMAX, including red leaks.
3452
C MAXK marks end of expected red-leak bands, starting at NBANDS+1.
3455
C EXAMPLE of red-leak storage:
3457
C In UBVRI system, NBANDS = 5. Suppose we measure only B,V,R, & BRL.
3462
C 2 B 6 <-- KMIN = 2
3464
C 4 R 0 <-- LASTK = 4
3466
C 6 BRL 0 <-- KMAX = MAXK = 6
3471
DO 1000 K=KMIN,LASTK
3472
IF (NBRL(K).EQ.0) GO TO 1000
3476
DO 900 NIGHT=1,NIGHTS
3479
C get limits for current night & passband.
3480
IF (NITE(J).NE.NIGHT) GO TO 20
3481
IF (DTYPE(J).EQ.'Y') GO TO 10
3489
C JBGN, JEND index the first & last data of ANY kind for the night.
3490
C JMIN, JMAX index the first & last STAR data for the night & band.
3491
C KK counts star data for the night.
3496
IF (DTYPE(J).EQ.'Y' .OR. NSTR(J).LT.0) GO TO 100
3497
C here if datum is stellar.
3499
IF (KBND(J).EQ.K) THEN
3501
C this observation requires red-leak correction.
3503
C scan forward & back for nearest red-leak datum for same star:
3507
IF (J1.GT.JMAX) GO TO 40
3508
IF (KBND(J1).NE.KRL .OR. NSTR(J1).NE.ISTR) GO TO 40
3516
IF (J2.LT.JMIN) GO TO 50
3517
IF (KBND(J2).NE.KRL .OR. NSTR(J2).NE.ISTR) GO TO 50
3523
CALL TV('No red-leak datum found for')
3524
CALL TVN(STARS(ISTR))
3529
CALL TV('Negative red-leak found for')
3530
CALL TVN(STARS(NSTR(J)))
3534
C Correct for red leak:
3535
SIGNAL(J)=SIGNAL(J)-RL*RLFACT(K)
3536
PHVAR(J)=PHVAR(J)+RLVAR*RLFACT(K)
3537
IF(ESTER(J).NE.0.)ESTER(J)=SQRT(ESTER(J)**2+ERVAR*RLFACT(K))
3543
C prepare for next pass.
3554
SUBROUTINE NEED(LINES)
3556
C Tell it how many lines you need for a plot. It reads the LOG
3557
C keyword and adds spaces to start a new page if necessary.
3562
INTEGER LINES, NGOT, LINE, IUNIT, NULLS, ISTAT, NLPP, I
3565
CALL STKRDI('LOG',6,1, NGOT,LINE,IUNIT,NULLS,ISTAT)
3566
IF (ISTAT.NE.0) CALL TERROR(9999,'Could not read Keyword LOG')
3567
CALL STKRDI('LOG',7,1, NGOT,NLPP,IUNIT,NULLS,ISTAT)
3568
IF (ISTAT.NE.0) CALL TERROR(9999,'Could not read Keyword LOG')
3569
IF(LINE+LINES.GT.NLPP)THEN
3570
C we must space to a new page.
added
removed
contrib/pepsys/libsrc/redsubs2.for
