C @(#)redsubs2.for 17.1.1.1 (ES0-DMD) 01/25/02 17:17:18 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 SUBROUTINE POLYGON(ARRAY,IFLAG,LABEL) C @(#)redsubs2.for 17.1.1.1 (ESO-IPG) 01/25/02 17:17:18 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT (c) European Southern Observatory 1992 C.IDENT redsubs2.for C.MODULE C.AUTHOR Andrew T. Young C.KEYWORD C.LANGUAGE FORTRAN 77 C.PURPOSE more subroutines for REDUCE C.COMMENTS C.VERSION 5.9 C.RETURNS C.ENVIRONMENT C. C----------------------------------------------------------------------------- C C C scans array for items with NSTR=iflag, and does polygon interpolation. C C IMPLICIT NONE REAL ARRAY, XBAR, YBAR, SLOPE, FRACT, YMAX, YMIN C INTEGER IFLAG, LEN, LWORD, JBGN, NIGHT, K, J, 1 JMIN, JMAX, JEND, K1, K2 C INCLUDE 'MID_REL_INCL:obs.inc' C C DIMENSION ARRAY(MXOBS) C CHARACTER A, LABEL*(*), CARD*79 C LOGICAL HELP EXTERNAL HELP C C C LEN=LWORD(LABEL) JBGN=1 C DO 100 NIGHT=1,NIGHTS C CALL SPACE2 CALL NEED(24) CARD(:LEN)=LABEL(:LEN) CARD(LEN+1:LEN+6)=' for ' CALL JD2DAT(REAL(2400000.5D0+DTZERO(NIGHT)),CARD(LEN+7:)) CALL TV(CARD) K=0 YMAX=-3.E33 YMIN=3.E33 C DO 10 J=JBGN,NDATA C C find data for current night. IF (NITE(J).NE.NIGHT) GO TO 20 IF (NSTR(J).EQ.IFLAG)THEN C add to XS, YS. K=K+1 XS(K)=DJOBS(J) YS(K)=ARRAY(J) YMAX=MAX(YMAX,YS(K)) YMIN=MIN(YMIN,YS(K)) IF (K.EQ.1) JMIN=J JMAX=J END IF C 10 CONTINUE J=NDATA+1 C 20 CONTINUE JEND=J-1 C C JBGN, JEND index the first and last data of ANY kind for the night. C JMIN, JMAX index the first and last IFLAG data for the night. C C C K is count of data for night N. C IF (YMAX.GT.YMIN) THEN C K > 1. 22 CALL NEED(24) CALL PLOT(K,XS,YS,'*') CALL TVN(' Time (decimal days) -->') CALL ASKN('OK to interpolate with a polygon?',A) IF (A.EQ.'N') THEN 24 CALL TV('Enter C to replace data with a Constant,') IF (K.GE.6) 1 CALL TVN(' L to smooth with a Linear fit,') CALL ASKN(' or Q to Quit.',A) IF (A.EQ.'C') THEN C be lazy: fit line, but force it horizontal. CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE) SLOPE=0. ELSE IF (A.EQ.'L') THEN IF (K.LT.6) THEN CALL TV('Not enough data to fit a line safely.') CALL TVN('Try a constant instead.') GO TO 24 END IF CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE) ELSE GO TO 24 END IF C DO 28 J=JBGN,JEND ARRAY(J)=(DJOBS(J)-XBAR)*SLOPE+YBAR 28 CONTINUE C GO TO 41 ELSE IF (A.EQ.'Y' .OR. A.EQ.'O') THEN C go on. ELSE IF (HELP(A)) THEN CALL TV('Reply NO to smooth data.') GO TO 22 ELSE GO TO 22 END IF ELSE IF (YMAX.EQ.YMIN) THEN CALL TV('Only 1 value available for this night!') CALL SPACE2 DO 30 J=JBGN,JEND ARRAY(J)=ARRAY(JMIN) 30 CONTINUE GO TO 99 ELSE CALL TV('NO DATA AVAILABLE FOR THIS NIGHT') CALL STSEPI END IF C C Now interpolate. C K1=1 K2=2 C DO 40 J=JBGN,JEND IF (J.LT.JMIN) THEN ARRAY(J)=ARRAY(JMIN) ELSE IF (J.GT.JMAX) THEN ARRAY(J)=ARRAY(JMAX) ELSE IF (DJOBS(J).GT.XS(K2)) THEN K1=K2 K2=K2+1 END IF C C linear interp... FRACT=(DJOBS(J)-XS(K1))/(XS(K2)-XS(K1)) ARRAY(J)=YS(K1)+(YS(K2)-YS(K1))*FRACT END IF 40 CONTINUE C 41 CARD(66:78)=CARD(53:65) CARD(53:65)=CARD(40:52) CARD(40:52)=CARD(27:39) CARD(27:39)=CARD(14:26) CARD(14:26)=CARD(:13) CARD(:13)='Interpolated ' CALL NEED(24) CALL TV(CARD) C fudge to get right plot limits.... XS(1)=DJOBS(JBGN) XS(K)=DJOBS(JEND) CALL PLOT(-K,XS,YS,' ') DO 50 J=JBGN,JEND CALL PLOT(-1,DJOBS(J),ARRAY(J),'+') 50 CONTINUE C now undo fudge. XS(1)=DJOBS(JMIN) XS(K)=DJOBS(JMAX) CALL PLOT(K,XS,YS,'*') CALL RTNCON(' Time (decimal days) -->',28) C 99 JBGN=JEND+1 C IF (JBGN.GT.NDATA) RETURN 100 CONTINUE C RETURN END SUBROUTINE DARKFIT(NDET,NDUSED,LABEL) C C scans SIGNAL for items with NSTR=-1, models dark level vs. DTMP, C and subtracts dark current for detector number NDET. C C C C A U T I O N : This routine has *** NOT *** been debugged! C =============== C IMPLICIT NONE C REAL XLIMS, YLIMS, SIGMIN, SIGMAX, YMAX, 1 YMIN, XBAR, YBAR, SLOPE, XBAR1, YBAR1, SLOPE1, TRY, XBAR2, 2 YBAR2, SLOPE2, FIT INTEGER NDET, LEN, LWORD, K, J, JMIN, 1 JMAX, K3, K1, I, K2, NIX C INCLUDE 'MID_REL_INCL:mbands.inc' C PARAMETER (MBANDS=9) INTEGER MXCOLR PARAMETER (MXCOLR=3*MBANDS) INTEGER NDUSED(MXCOLR) C INCLUDE 'MID_REL_INCL:obs.inc' C This defines SIGNAL, NSTR, KBND, etc. C C Includes for dlsq: INCLUDE 'MID_REL_INCL:kingfit.inc' INCLUDE 'MID_REL_INCL:dlsq.inc' C CHARACTER A1, LABEL*16, CARD*79 C DIMENSION XLIMS(2),YLIMS(2) C LOGICAL HELP EXTERNAL HELP C EXTERNAL YP2EXP C C C Begin execution: C LEN=LWORD(LABEL) C 2 K=0 SIGMAX=-3.E33 SIGMIN=3.E33 CALL SPACE2 CARD(:LEN)=LABEL(:LEN) CARD(LEN+1:)=' Log(DARK) vs. Temperature' C DO 10 J=1,NDATA C C find dark data for current detector. IF (NSTR(J).EQ.-1 .AND. KBND(J).EQ.NDET)THEN C add to XS, YS. K=K+1 XS(K)=DTMP(J) IF (SIGNAL(J).GT.0.) THEN YS(K)=LOG(SIGNAL(J)) ELSE CARD='Negative dark reading at MJD =' WRITE(CARD(35:),'(F12.5)') DTZERO(NITE(J))+DJOBS(J) CARD(48:)='in channel' WRITE(CARD(55:),'(I3)') NDET CALL TV(CARD) K=K-1 GO TO 10 END IF SIGMAX=MAX(SIGMAX,SIGNAL(J)) SIGMIN=MIN(SIGMIN,SIGNAL(J)) C C SYMB is night label: IF (NITE(J).LT.10) THEN SYMB(K)=CHAR(ICHAR('0')+NITE(J)) ELSE SYMB(K)=CHAR(ICHAR('A')+NITE(J)-1) END IF C IF (K.EQ.1) THEN JMIN=J END IF C do this every time, in case only 1 datum. JMAX=J END IF C 10 CONTINUE C C K is count of data for current detector. C 21 YMAX=-3.E33 YMIN=3.E33 C IF (SIGMAX.GT.SIGMIN) THEN C K > 1. C set page size: CALL PLOT(0,78.,22.,'P') C set different-symbol mode: 22 CALL PLOT(0,78.,22.,'D') CALL NEED(24) CALL TV(CARD) C plot each night with its own symbol. CALL PLOT(K,XS,YS,SYMB) C set same-symbol mode: CALL PLOT(0,78.,22.,'S') IF (K.GT.10) THEN CALL TVN(' Temperature -->') CALL ASKN('OK to fit with 2 lines?',A1) ELSE IF (HELP(A1)) THEN CALL TV('Reply YES only if concave upward') CALL TVN('(high at ends, lower in middle)') GO TO 22 ELSE A1='N' CALL RTNCON(' Temperature -->', 20) END IF C IF (A1.EQ.'N') THEN 24 CALL TV('Enter C to replace data with a Constant,') CALL TVN(' L to smooth with a Linear fit,') CALL ASKN(' or Q to Quit.',A1) C IF (A1.EQ.'C') THEN C be lazy: fit line, but force it horizontal. CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE) SLOPE=0. NIX=1 IX(1)=2 ELSE IF (A1.EQ.'L') THEN IF (K.LT.6) THEN CALL TV('Not enough data to fit 1 line safely.') CALL TVN('Try a constant instead.') GO TO 24 END IF C fit a simple exponential in temperature: CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE) NIX=0 ELSE GO TO 24 END IF C C C Refine one-exponential fit: C IPR=0 KIP=2 PG(1)=EXP(YBAR) PG(2)=SLOPE DO 26 J=1,K X(1,J)=XS(J)-XBAR Y(J)=EXP(YS(J)) C use expected ordinates to weight: W(J)=1./EXP(YBAR+X(1,J)*SLOPE) 26 CONTINUE C refine fit: C subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST) CALL DLSQ(YP2EXP,K,2, 1, 2, NIX, 1, 0, 1.E-3) SLOPE=LOG(P(2)) YBAR=P(1) C C C now subtract fit from ALL data with this detector: C DO 28 J=1,NDATA IF ((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR. 1 (NSTR(J).EQ.-1 .AND. KBND(J).NE.NDET)) GO TO 28 C subtract dark: SIGNAL(J)=SIGNAL(J)-EXP((DTMP(J)-XBAR)*SLOPE+YBAR) IF (NSTR(J).EQ.-1) THEN YMAX=MAX(YMAX,SIGNAL(J)) YMIN=MIN(YMIN,SIGNAL(J)) END IF 28 CONTINUE C GO TO 61 C C ELSE IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN C go on to do two-exponential fit at *****. C ELSE IF (A1.EQ.'H')THEN CALL TV('If the plot is generally concave upward,') CALL TVN('and the different nights (different symbols)') CALL TVN('appear to follow the same curve, reply YES.') CALL TV('If the plot is straight or concave downward,') CALL TVN('or some nights are offset, reply NO.') CALL TV('Overlapping points are denoted by "$".') CALL SPACE2 CALL RTNCON(' ',1) GO TO 21 C ELSE GO TO 22 C END IF C ELSE IF (SIGMAX.EQ.SIGMIN) THEN CALL TV('Only 1 dark value available!') CALL SPACE2 DO 30 J=1,NDATA C subtract dark: IF (NSTR(J).GT.0 .AND. NDUSED(KBND(J)).EQ.NDET) 1 SIGNAL(J)=SIGNAL(J)-SIGMIN 30 CONTINUE C skip plot. GO TO 99 ELSE CALL TV('NO DATA AVAILABLE -- cannot subtract dark') RETURN END IF C C ***** C C Here to fit 2 lines (exponentials), and remove their sum. C K3=K/3 C IF (K3.LT.9) THEN CALL TV('CAUTION: probably too few points to get a fit.') CALL ASKN('Were there any wild points in the plot above?',A1) C IF (A1.EQ.'Y')THEN CALL TV('Then you should try to fit just one line.') CALL RTNCON(' ',1) GO TO 21 ELSE IF (A1.EQ.'N') THEN C go ahead. ELSE IF (A1.EQ.'H') THEN C get help: show plot again. GO TO 21 ELSE END IF C END IF C CALL SORT2(XS,YS,K) C fit trial line to bottom 1/3 of data. CALL ROBLIN(XS,YS,K3, XBAR1,YBAR1,SLOPE1) C C remove lower trial line, and fit the upper line. C K1=0 C DO 32 I=1,K3 C move upper third to lower third: K2=K-K3+I TRY=EXP(YS(K2))-EXP(YBAR1+SLOPE1*(XS(K2)-XBAR1)) IF (TRY.LE.0.) THEN C skip this point. ELSE K1=K1+1 C save K1 in K1+K. XS(K1+K)=XS(K1) YS(K1+K)=YS(K1) C move K2 to K1. YS(K1)=LOG(TRY) XS(K1)=XS(K2) END IF 32 CONTINUE C IF (K1.LT.6)THEN CALL TV('Sorry; can''t fit 2 lines. Try 1 line.') CALL RTNCON(' ',1) GO TO 2 ELSE CALL ROBLIN(XS,YS,K1, XBAR2,YBAR2,SLOPE2) IF (SLOPE2.LT.0. .OR. SLOPE2.LT.SLOPE1) THEN C unphysical result; complain. CALL TV('Sorry; attempted fit is nonsense.') CALL RTNCON('Try again, using 1 line.',24) GO TO 2 ELSE END IF END IF C C Now restore saved data to lower K1: DO 35 I=1,K1 XS(K1)=XS(K1+K) YS(K1)=YS(K1+K) 35 CONTINUE C C ***** Polish the result. ***** C C Set up for DLSQ: C IPR=0 NIX=0 C Set up data: XBAR=(XBAR1+XBAR2)*0.5 DO 36 I=1,K X(1,I)=XS(I)-XBAR Y(I)=YS(I) W(I)=1./(EXP(YBAR1+SLOPE1*(XS(I)-XBAR1)) + 1 EXP(YBAR2+SLOPE2*(XS(I)-XBAR2)) ) 36 CONTINUE KIP=4 PG(1)=EXP(YBAR1+SLOPE1*(XBAR-XBAR1)) PG(2)=SLOPE1 PG(3)=EXP(YBAR2+SLOPE2*(XBAR-XBAR2)) PG(4)=SLOPE2 C refine fit: DO 40 I=1,3 C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST) CALL DLSQ(YP2EXP, K, 2, 1, 2, NIX, 1, 0, 1.E-3) C IF (P(1).LE.0.D0 .OR. P(3).LE.0.D0) THEN CALL TV('Sorry -- unable to fit two lines. Try one.') CALL RTNCON(' ',1) GO TO 2 END IF C Update values: DO 38 J=1,4 38 PG(J)=P(J) SLOPE1=P(2) YBAR1=LOG(P(1))-SLOPE1*(XBAR-XBAR1) SLOPE2=P(4) YBAR2=LOG(P(3))-SLOPE2*(XBAR-XBAR2) C Revise weights: DO 39 J=1,K W(J)=1./(EXP(YBAR1+SLOPE1*(XS(I)-XBAR1)) + 1 EXP(YBAR2+SLOPE2*(XS(I)-XBAR2)) ) 39 CONTINUE 40 CONTINUE C C C C Subtract the interpolated values from ALL data for this det. C DO 60 J=1,NDATA IF ((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR. 1 (NSTR(J).EQ.-1 .AND. KBND(J).NE.NDET)) GO TO 60 FIT=EXP(YBAR1+SLOPE1*(DTMP(J)-XBAR1)) + 1 EXP(YBAR2+SLOPE2*(DTMP(J)-XBAR2)) SIGNAL(J)=SIGNAL(J)-FIT IF (NSTR(J).NE.-1) GO TO 60 C use only dark data to set limits for residual plot. YMAX=MAX(YMAX,SIGNAL(J)) YMIN=MIN(YMIN,SIGNAL(J)) 60 CONTINUE C C Now inspect resids. vs. time: C 61 CARD(:LEN)=LABEL(:LEN) CARD(LEN+1:)=' Dark RESIDUALS vs. Time' CALL NEED(24) CALL TV(CARD) C set plot limits.... XLIMS(1)=DJOBS(1) XLIMS(2)=DJOBS(NDATA) YLIMS(1)=YMIN YLIMS(2)=YMAX CALL PLOT(0,XLIMS,YLIMS,'L') DO 80 J=1,NDATA IF (NSTR(J).NE.-1 .OR. KBND(J).NE.NDET) GO TO 80 IF (NITE(J).LT.10) THEN A1=CHAR(ICHAR('0')+NITE(J)) ELSE A1=CHAR(ICHAR('A')+NITE(J)-1) END IF CALL PLOT(-1,DJOBS(J),SIGNAL(J),A1) 80 CONTINUE C force plot: CALL PLOT(1,0.,0.,'+') CALL RTNCON(' Time (decimal days) -->',28) C 99 CONTINUE C RETURN END SUBROUTINE YP2EXP C C Fits Y = P(1)*EXP(P(2)*Z(1)) + P(3)*EXP(P(4)*Z(1)) C C IMPLICIT NONE C C Includes for dlsq: INCLUDE 'MID_REL_INCL:kingfit.inc' INCLUDE 'MID_REL_INCL:dlsq.inc' C DOUBLE PRECISION TERM C C C Start with first term: C C C KIP is set by DARKFIT. C IP(1)=1 IP(2)=2 C TERM=P(2)*Z(1) TERM=EXP(MIN(TERM,30.D0)) PART(1)=TERM PART(2)=P(1)*TERM*Z(1) C YT = P(1)*TERM C IF (KIP.EQ.2) RETURN C C Here for second term: C IP(3)=3 IP(4)=4 C TERM=P(4)*Z(1) TERM=EXP(MIN(TERM,30.D0)) PART(3)=TERM PART(4)=P(3)*TERM*Z(1) C YT = YT + P(3)*TERM C C RETURN END SUBROUTINE DARKPOLY(NDET,NDUSED,LABEL) C C scans SIGNAL for items with NSTR=-1, and does polygon interpolation C and subtraction of dark current for detector number NDET. C C IMPLICIT NONE C REAL XLIMS, YLIMS, SIGMIN, SIGMAX, SIGA, SIGZ, YMAX, 1 YMIN, XBAR, YBAR, SLOPE, FRACT INTEGER NDET, LEN, LWORD, JBGN, NIGHT, K, J, 1 JMIN, JMAX, JEND, K1, K2, NIX C INCLUDE 'MID_REL_INCL:mbands.inc' C PARAMETER (MBANDS=9) INTEGER MXCOLR PARAMETER (MXCOLR=3*MBANDS) INTEGER NDUSED(MXCOLR) C INCLUDE 'MID_REL_INCL:obs.inc' C This defines NITE, SIGNAL, NSTR, KBND, etc. C C Includes for dlsq: INCLUDE 'MID_REL_INCL:kingfit.inc' INCLUDE 'MID_REL_INCL:dlsq.inc' C CHARACTER A1, LABEL*16, CARD*79 C DIMENSION XLIMS(2),YLIMS(2) C LOGICAL HELP EXTERNAL HELP C EXTERNAL YPLIN C C C Begin execution: C LEN=LWORD(LABEL) JBGN=1 C DO 100 NIGHT=1,NIGHTS C CALL SPACE2 CALL NEED(24) CARD(:LEN)=LABEL(:LEN) CARD(LEN+1:LEN+10)=' DARK on ' CALL JD2DAT(REAL(2400000.5D0+DTZERO(NIGHT)),CARD(LEN+11:)) CARD(LWORD(CARD)+3:)='= MJD' WRITE (CARD(LWORD(CARD)+1:),'(F7.0)') DTZERO(NIGHT) CALL TV(CARD) SIGMAX=-3.E33 SIGMIN=3.E33 K=0 C DO 10 J=JBGN,NDATA C C find dark data for current night and detector. IF (NITE(J).NE.NIGHT) GO TO 20 c KBND holds NDET for DARK data, remember. IF (NSTR(J).EQ.-1 .AND. KBND(J).EQ.NDET)THEN C add to XS, YS. K=K+1 XS(K)=DJOBS(J) YS(K)=SIGNAL(J) IF (K.EQ.1) JMIN=J IF (K.EQ.1) SIGA=SIGNAL(J) JMAX=J SIGZ=SIGNAL(J) SIGMAX=MAX(SIGMAX,SIGNAL(J)) SIGMIN=MIN(SIGMIN,SIGNAL(J)) END IF C 10 CONTINUE J=NDATA+1 C 20 CONTINUE JEND=J-1 C C JBGN, JEND index the first and last data of ANY kind for the night. C JMIN, JMAX index the first and last DARK data for the night. C SIGA, SIGZ are the first and last dark LEVELS for the night. C SIGMIN, SIGMAX are the lowest and highest dark LEVELS for the night. C C YMAX=-3.E33 YMIN=3.E33 C C K is count of data for night N. C IF (SIGMAX.GT.SIGMIN) THEN 22 CALL NEED(24) CALL PLOT(K,XS,YS,'*') CALL TVN(' TIME (decimal days) -->') CALL ASKN('OK to interpolate with a polygon?',A1) IF (A1.EQ.'N') THEN 24 CALL TV('Enter C to replace data with a Constant,') IF (K.GE.6) 1 CALL TVN(' L to smooth with a Linear fit,') CALL ASKN(' or Q to Quit.',A1) IF (A1.EQ.'C') THEN C be lazy: fit line, but force it horizontal. CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE) SLOPE=0. NIX=1 IX(1)=2 ELSE IF (A1.EQ.'L') THEN IF (K.LT.6) THEN CALL TV('Not enough data to fit a line safely.') CALL TVN('Try a constant instead.') GO TO 24 END IF CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE) NIX=0 ELSE IF (HELP(A1)) THEN CALL TV('Reply NO to smooth data.') GO TO 22 ELSE GO TO 24 END IF C C C Refine fit: C IPR=0 PG(1)=YBAR PG(2)=SLOPE DO 26 J=1,K X(1,J)=XS(J)-XBAR Y(J)=YS(J) W(J)=1. 26 CONTINUE C refine fit: C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST) CALL DLSQ(YPLIN,K,2, 1, 2, NIX, 1, 0, 1.E-3) SLOPE=P(2) YBAR=P(1) C C C Subtract fit: C DO 28 J=JBGN,JEND IF((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR. C real star, not observed with current det... 1 (NSTR(J).LE.0 .AND. KBND(J).NE.NDET)) GO TO 28 C ...or dark, not observed with current det. C OK, subtract interpolated dark line: SIGNAL(J)=SIGNAL(J)-((DJOBS(J)-XBAR)*SLOPE+YBAR) IF (NSTR(J).EQ.-1) THEN C reset limits for residual plot. YMAX=MAX(YMAX,SIGNAL(J)) YMIN=MIN(YMIN,SIGNAL(J)) END IF 28 CONTINUE C C plot residuals. C CALL NEED(24) CARD(66:78)=CARD(53:65) CARD(53:65)=CARD(40:52) CARD(40:52)=CARD(27:39) CARD(27:39)=CARD(14:26) CARD(14:26)=CARD(:13) CARD(:13)='RESIDUALS of ' CALL TV(CARD) C set plot limits.... XLIMS(1)=DJOBS(JBGN) XLIMS(2)=DJOBS(JEND) YLIMS(1)=YMIN YLIMS(2)=YMAX CALL PLOT(0,XLIMS,YLIMS,'L') CALL XAXIS(XLIMS) DO 29 J=JBGN,JEND IF (NSTR(J).NE.-1 .OR. KBND(J).NE.NDET) GO TO 29 C plot only dark resids., not star data. CALL PLOT(-1,DJOBS(J),SIGNAL(J),'*') 29 CONTINUE C force plot: CALL PLOT(1,0.,0.,'+') CALL RTNCON(' Time (decimal days) -->',28) GO TO 99 C ELSE IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN C go on. ELSE GO TO 22 END IF ELSE IF (K.EQ.1) THEN CALL TV('Only 1 datum available for this night!') CALL SPACE2 DO 30 J=JBGN,JEND C subtract dark: IF((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR. 1 (NSTR(J).LE.0 .AND. KBND(J).NE.NDET)) THEN c skip it. ELSE SIGNAL(J)=SIGNAL(J)-SIGA END IF 30 CONTINUE C skip plot. GO TO 99 ELSE CALLTV('NO DATA AVAILABLE FOR THIS NIGHT -- no subtraction') GO TO 99 END IF C C Do polygon interpolation here. C K1=1 K2=2 C DO 40 J=JBGN,JEND IF (NDUSED(KBND(J)).NE.NDET) GO TO 40 C skip observations with other detectors. IF (J.LT.JMIN) THEN C Before 1st dark (at JMIN); use it. SIGNAL(J)=SIGNAL(J)-SIGA ELSE IF (J.GT.JMAX) THEN C After last dark (at JMAX); use it. SIGNAL(J)=SIGNAL(J)-SIGZ ELSE C Can interpolate. IF (DJOBS(J).GT.XS(K2)) THEN C Step to next interval. K1=K2 K2=K2+1 END IF C C linear interp... FRACT=(DJOBS(J)-XS(K1))/(XS(K2)-XS(K1)) SIGNAL(J)=SIGNAL(J)-(YS(K1)+(YS(K2)-YS(K1))*FRACT) END IF C note that we have no plot of residuals at the DARK data points, C as they are negligible for polygon fit. 40 CONTINUE C 99 JBGN=JEND+1 C IF (JBGN.GT.NDATA) RETURN 100 CONTINUE C RETURN END SUBROUTINE ROBLIN (X, Y, N, XBAR, YBAR, SLOPE) C C Fits a robust line to N pairs (X,Y), having medians (Xbar, Ybar) C and SLOPE. Method: see Hoaglin, Mosteller & Tukey, Ch.5. C C Does NOT mess up input data; but leaves them sorted on X in XDUM,YDUM. C C C Input arguments: C INTEGER N REAL X(N), Y(N) C C Output arguments: C REAL XBAR,YBAR,SLOPE C C Scratch arrays: C INCLUDE 'MID_REL_INCL:obs.inc' C (defines MXOBS.) REAL XDUM,YDUM,ZDUM COMMON /SCRAT/ XDUM(MXOBS), YDUM(MXOBS), ZDUM(MXOBS) C C Local variables: C INTEGER I,J,N3,MID,MIDL,MIDR REAL XL,XR,YL,YR,B,B0,B1,B2,RL,RR,DELTA1,DELTAX,DR0,DR1 C C C ***** BEGIN EXECUTION ***** C IF (N.GT.MXOBS)THEN CALL TV('Too many data to sort in available slots.') CALL TV('Expand MXOBS and recompile.') CALL STSEPI ELSE IF (N.EQ.2)THEN XBAR=(X(1)+X(2))/2. YBAR=(Y(1)+Y(2))/2. IF (X(2).NE.X(1)) THEN SLOPE=(Y(2)-Y(1))/(X(2)-X(1)) ELSE SLOPE=0. END IF RETURN ELSE IF (N.EQ.1)THEN XBAR=X(1) YBAR=Y(1) SLOPE=0. RETURN ELSE IF (N.LE.0)THEN CALL TV('ROBLIN can''t fit a line to 0 points!') CALL TVN('BUG in program') CALL STSEPI END IF C DO 1 I=1,N XDUM(I)=X(I) YDUM(I)=Y(I) ZDUM(I)=Y(I) 1 CONTINUE C CALL SORT1(ZDUM,N) CALL SORT2(XDUM,YDUM,N) MID=(N+1)/2 C IF (MID.EQ.N/2) THEN C N EVEN XBAR=(XDUM(MID)+XDUM(MID+1))*0.5 YBAR=(ZDUM(MID)+ZDUM(MID+1))*0.5 ELSE C N ODD XBAR=XDUM(MID) YBAR=ZDUM(MID) END IF C C Pick size of end thirds: C N3=(N+1)/3 C Pick end thirds: MIDL=(N3+1)/2 MIDR=N+1-MIDL C IF (MIDL.EQ.N3/2) THEN C N3 EVEN XL=(XDUM(MIDL)+XDUM(MIDL+1))*0.5 XR=(XDUM(MIDR)+XDUM(MIDR-1))*0.5 ELSE C N3 ODD XL=XDUM(MIDL) XR=XDUM(MIDR) END IF DELTAX=(XR-XL) C C Now find y medians for left and right groups. C LEFT group: DO 10 I=1,N3 10 ZDUM(I)=YDUM(I) CALL SORT1(ZDUM,N3) C IF (MIDL.EQ.N3/2) THEN C N3 EVEN YL=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5 ELSE C N3 ODD YL=ZDUM(MIDL) END IF C RIGHT group: DO 20 I=1,N3 20 ZDUM(I)=YDUM(N+1-I) CALL SORT1(ZDUM,N3) C IF (MIDL.EQ.N3/2) THEN C N3 EVEN YR=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5 ELSE C N3 ODD YR=ZDUM(MIDL) END IF C B0=(YR-YL)/DELTAX B=B0 C C B0 is initial slope. C C ***** BEGIN ITERATIONS ***** DO 50 J=1,5 C Put RESIDS. in zdum: C LEFT: DO 30 I=1,N3 30 ZDUM(I)=YDUM(I)-(YBAR+B*(XDUM(I)-XBAR)) CALL SORT1(ZDUM,N3) C IF (MIDL.EQ.N3/2) THEN C N3 EVEN RL=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5 ELSE C N3 ODD RL=ZDUM(MIDL) END IF C RIGHT: DO 40 I=1,N3 40 ZDUM(I)=YDUM(N+1-I)-(YBAR+B*(XDUM(N+1-I)-XBAR)) CALL SORT1(ZDUM,N3) C IF (MIDL.EQ.N3/2) THEN C N3 EVEN RR=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5 ELSE C N3 ODD RR=ZDUM(MIDL) END IF C IF (J.EQ.1) THEN DR0=RR-RL DELTA1=DR0/DELTAX B1=B0+DELTA1 IF (B1.EQ.B0) GOTO 100 B=B1 ELSE DR1=RR-RL IF (DR1.EQ.0.) GOTO 100 IF (ABS((DR1-DR0)/DR1).LT.1.E-4) GOTO 100 B2=B1-DR1*(B1-B0)/(DR1-DR0) B=B2 IF (ABS(B2-B1).LT.1.E-6*ABS(B1)) GO TO 100 C prepare for next round: B0=B1 B1=B2 DR0=DR1 END IF 50 CONTINUE C 100 CONTINUE SLOPE=B C Current value of b is converged. Recalculate all resids.: C DO 110 I=1,N 110 ZDUM(I)=YDUM(I)-(YBAR+B*(XDUM(I)-XBAR)) CALL SORT1(ZDUM,N) C IF (MID.EQ.N/2) THEN C N EVEN YBAR=(ZDUM(MID)+ZDUM(MID+1))*0.5+YBAR ELSE C N ODD YBAR=ZDUM(MID)+YBAR END IF C C RETURN END SUBROUTINE SORT1(X,N) C C Sorts X, in place. C C BOOTHROYD'S SHELL-SORT: COMM.ACM 6,445(1963) ALG.201. C C SEE H.LORIN,"SORTING & SORT SYSTEMS" (ADDISON-WESLEY,1975) C CH.3 & 11 FOR THEORY AND TESTS. C C IMPLICIT NONE C INTEGER N REAL X(N) C C C Local variables: C INTEGER I,J,M REAL DUM C C **** BEGIN **** C IF(N.EQ.1)RETURN I=1 3 I=I+I IF(I.LE.N)GO TO 3 M=I-1 5 M=M/2 IF(M.EQ.0)RETURN DO 20 J=1,N-M DO 10 I=J,1,-M IF(X(I+M).GE.X(I)) GO TO 20 DUM=X(I) X(I)=X(I+M) X(I+M)=DUM 10 CONTINUE 20 CONTINUE GO TO 5 C END SUBROUTINE SORT2(X,Y,N) C C Sorts (X,Y) pairs on X, in place. C C BOOTHROYD'S SHELL-SORT: COMM.ACM 6,445(1963) ALG.201. C C SEE H.LORIN,"SORTING & SORT SYSTEMS" (ADDISON-WESLEY,1975) C CH.3 & 11 FOR THEORY AND TESTS. C C IMPLICIT NONE C INTEGER N REAL X(N),Y(N) C C C Local variables: C INTEGER I,J,M REAL DUM C C **** BEGIN **** C IF(N.EQ.1)RETURN I=1 3 I=I+I IF(I.LE.N)GO TO 3 M=I-1 5 M=M/2 IF(M.EQ.0)RETURN DO 20 J=1,N-M DO 10 I=J,1,-M IF(X(I+M).GE.X(I)) GO TO 20 DUM=X(I) X(I)=X(I+M) X(I+M)=DUM DUM=Y(I) Y(I)=Y(I+M) Y(I+M)=DUM 10 CONTINUE 20 CONTINUE GO TO 5 C END SUBROUTINE YPLIN C C YP for use by DLSQ in linear fits. C C C fits: Y = P(1) + P(2)*Z(1) C C IMPLICIT NONE C C Includes for dlsq: INCLUDE 'MID_REL_INCL:kingfit.inc' INCLUDE 'MID_REL_INCL:dlsq.inc' C C C YT = P(1) + P(2)*Z(1) C PART(1)=1.D0 PART(2)=Z(1) C IP(1)=1 IP(2)=2 KIP=2 C C RETURN END SUBROUTINE SKYSUB(BANDS,KMIN,KMAX,DTYPE,KDIAPHRAGM,EXTIN) C C scans SIGNAL for items with DTYPE='Y', and does modelling, C interpolation, and subtraction of sky for each passband. C C IMPLICIT NONE C C INCLUDE 'MID_REL_INCL:mbands.inc' C PARAMETER (MBANDS=9) C INCLUDE 'MID_REL_INCL:obs.inc' REAL XDUM,YDUM,ZDUM COMMON /SCRAT/ XDUM(MXOBS), YDUM(MXOBS), ZDUM(MXOBS) C C C Args: C INTEGER KMIN, KMAX, KDIAPHRAGM CHARACTER DTYPE(MXOBS) REAL EXTIN(MBANDS) C C REAL COSOB, SINOB, RASUN, DESUN, HASUN, SOLONG, ELMOON, BMOON, 1 SINPHI, COSPHI, ELONG, ELROT, UTROT, STUTZ, TWOPI, ST2UT, 2 TNOON1, TNOON2, PI, DEGRAD, ALAT, XLIMS, YLIMS, ALTS, T, 3 YMAX, YMIN, AIRMAX, SALTS, OBLIQ, SIGMIN, 4 SIGMAX, RISE, SET, DAWN, DUSK, UTRAD, TLST, HAMOON, COSDEL, 5 SINDEL, COSHA, XM, YM, ZM, AZMOON, DELAZ, COSEL, XBAR, 6 YBAR, SLOPE, OLDVAR, ELONGMAX, ELONGMIN, ZENITH C INTEGER NB, LEN, LWORD, JBGN, NIGHT, NTWI, 1 K, J, JMIN, JMAX, JEND, MOONUP, NERROR, NIX, NPTS, NFIT C C COMMONS FOR SPHERICAL TRIG.: C COMMON /CSUN/ COSOB,SINOB,RASUN,DESUN,HASUN,SOLONG,ELMOON,BMOON COMMON /SPHERE/ SINPHI,COSPHI,ELONG,ELROT,UTROT,STUTZ,TWOPI,ST2UT, 1 TNOON1,TNOON2,PI,DEGRAD,ALAT INTEGER MXCOLR PARAMETER (MXCOLR=3*MBANDS) C C commons for star catalog: INCLUDE 'MID_REL_INCL:mstars.inc' C PARAMETER (MSTARS=1650) CHARACTER *32 STARS COMMON /SCATA/ STARS(MSTARS) C C Includes for dlsq: INCLUDE 'MID_REL_INCL:kingfit.inc' INCLUDE 'MID_REL_INCL:dlsq.inc' C C C Commons for 12 Diaphragms: C INTEGER NP2NDIA(12*MBANDS) INTEGER NDIA2NP(12,MBANDS), NDIAS,NDIA C indices: NDIA,K REAL XNDIA(12,MBANDS) COMMON /DIACOM/XNDIA,NP2NDIA,NDIA2NP,NDIAS,NDIA C CHARACTER DIANAM(12)*4 COMMON /DIACH/DIANAM C C C LOCAL variables: C CHARACTER A1, CARD*80, BANDS(MXCOLR)*8, DATSTR*16 C REAL ERRSAVE(12) INTEGER LORDER(12),L,NN, JDUSK,JDAWN, NPARMS,NREFDIA C LOGICAL TWILIT, SHORT C C EXTERNAL YPLIN, YPSKY C C DIMENSION XLIMS(2),YLIMS(2) DIMENSION ALTS(2), SALTS(5) C C true altitude of Moon at rise & set: DATA ALTS/-.0145,-.0145/ C C C true altitude of Sun at limit of astronomical twilight: DATA SALTS/-.309,0.,0.,0.,-.309/ C C C C Begin execution: C CALL SPACE2 CALL TV(' SKY subtraction:') CALL TV('You will see graphs of sky brightness, plotted against') CALL TVN('time, airmass, and lunar elongation, for each band,') CALL TVN('on each night separately.') CALL TV('Then you can choose the method of fitting or') CALL TVN('interpolating sky brightness, for each case.') CALL SPACE CALL RTNCON(' ',1) CALL PLOT(0,78.,23.,'P') C DO 1000 NB=KMIN,KMAX LEN=LWORD(BANDS(NB)) JBGN=1 DO 1000 NIGHT=1,NIGHTS CALL SPACE2 CALL NEED(24) 2 CARD(:LEN)=BANDS(NB) CARD(LEN+1:LEN+10)=' SKY on' CALL JD2DAT(REAL(2400000.5D0+DTZERO(NIGHT)),CARD(LEN+11:)) DATSTR=CARD(LEN+6:) SHORT=.FALSE. C 3 K=0 DO 5 NDIA=1,NDIAS NDIA2NP(NDIA,1)=0 5 CONTINUE YMAX=-3.E33 YMIN=3.E33 AIRMAX=1. C DO 10 J=JBGN,NDATA C C find range of sky data for current night and passband. IF (NITE(J).NE.NIGHT) GO TO 12 IF (KBND(J).NE.NB) GO TO 10 IF (DTYPE(J).EQ.'Y')THEN C add to XS, YS. K=K+1 XS(K)=DJOBS(J) YS(K)=SIGNAL(J) YMAX=MAX(YMAX,SIGNAL(J)) YMIN=MIN(YMIN,SIGNAL(J)) AIRMAX=MAX(AIRMAX,AIRM(J)) IF (K.EQ.1) JMIN=J IF (K.EQ.1) SIGMIN=SIGNAL(J) JMAX=J SIGMAX=SIGNAL(J) END IF C 10 CONTINUE J=NDATA+1 C 12 CONTINUE JEND=J-1 C C JBGN, JEND index the first and last data of ANY kind for the night. C JMIN, JMAX index the first and last SKY data for the night & band. C SIGMIN, SIGMAX are the first and last sky LEVELS for the night. C C K is count of data for night N. C T=(DTZERO(NIGHT)+INT(DJOBS((JBGN+JEND)/2))-51544.5)/36525. CALL STUTZR(T) C STUTZR sets sidereal time for later use by UTMOON & UTSUN. C Set OBLIQ for SUN (called by UTSUN). OBLIQ=(23.43929-T*(1.30042E-2+T*(1.639E-7-T*5.036E-7)))*DEGRAD COSOB=COS(OBLIQ) SINOB=SIN(OBLIQ) C C Adjust YMAX for possible outliers: C CALL SORT2(YS,XS,K) C (Note reversal of Y and X.) YM=(YS((K+1)/2)+YS(K/2+1))*0.5 C YM is median. NN=(K+2)/4 ZM=YS(K-NN)-YS(NN) C ZM is interquartile range. YMAX=MIN(YMAX,YM+4.*ZM) C IF (K.GT.1) THEN C C Take shortcut if required: IF (SHORT) THEN A1='T' GO TO 102 END IF C C set up plots; first, TIME: C C set plot limits: XLIMS(1)=DJOBS(JMIN) XLIMS(2)=DJOBS(JMAX) YLIMS(1)=YMIN YLIMS(2)=YMAX CALL PLOT(0,XLIMS,YLIMS,'L') C C plot moonrise & moonset. First, guarantee RISE < SET: RISE=DJOBS(JMIN) SET=DJOBS(JMAX) C UTMOON needs T in Julian centuries from J2000. T=((2400000.5D0-2451545.D0)+DTZERO(NIGHT)+ 1 INT(DJOBS((JBGN+JEND)/2)))/36525. CALL UTMOON(T,ALTS,1,*21) C UTMOON returns result in UTROT, in /SPHERE/. RISE=UTROT+INT(DJOBS((JBGN+JEND)/2)) C we have at most 18 intervals between YMIN & YMAX. CALL PLOT(-1,RISE,(YMIN+14.*YMAX)/15.,'|') CALL PLOT(-1,RISE,(YMIN*2.+13.*YMAX)/15.,'|') CALL PLOT(-1,RISE,YMAX,'R') IF (DJOBS(JMIN).LT.RISE .AND. DJOBS(JMAX).GT.RISE) 1 CARD(LWORD(CARD)+5:)='R = moonrise' 21 CALL UTMOON(T,ALTS,2,*22) SET=UTROT+INT(DJOBS((JBGN+JEND)/2)) CALL PLOT(-1,SET,(YMIN+14.*YMAX)/15.,'|') CALL PLOT(-1,SET,(YMIN*2.+13.*YMAX)/15.,'|') CALL PLOT(-1,SET,YMAX,'S') IF (DJOBS(JMIN).LT.SET .AND. DJOBS(JMAX).GT.SET) 1 CARD(LWORD(CARD)+5:)='S = moonset' 22 CONTINUE C C RISE & SET now are correct for Moon. Do SUN: DUSK=DJOBS(JMIN) DAWN=DJOBS(JMAX) C Here the argument is DAYS from J2000. CALL UTSUN(T*36525.,SALTS,1,*23) DUSK=UTROT+INT(DJOBS((JBGN+JEND)/2)) 23 CALL UTSUN(T*36525.,SALTS,5,*24) DAWN=UTROT+INT(DJOBS((JBGN+JEND)/2)) 24 CONTINUE C DUSK & DAWN are now limits of astronomical twilight. C TWILIT=.FALSE. JDUSK=JEND JDAWN=JBGN DO 26 J=JMIN,JMAX IF (KBND(J).NE.NB) GO TO 26 IF ( (DAWN.LT.DUSK .AND. 1 DJOBS(J).GT.DAWN .AND. DJOBS(J).LT.DUSK) .OR. 2 (DAWN .GT.DUSK .AND. 3 (DJOBS(J).GT.DAWN .OR. DJOBS(J).LT.DUSK) ) ) THEN C Twilight. A1='t' TWILIT=.TRUE. ELSE C Not twilit. A1='*' JDUSK=MIN(JDUSK,J) JDAWN=MAX(JDAWN,J) END IF IF (DTYPE(J).EQ.'Y')THEN CALL PLOT(-1,DJOBS(J),SIGNAL(J),A1) END IF 26 CONTINUE IF (TWILIT) CARD(LWORD(CARD)+5:)='t = twilight sky' CALL TV(CARD) CALL SPACE CALL PLOT(1,0.,0.,'+') WRITE(CARD,'(A,F7.0)')' TIME (MJD =',DTZERO(NIGHT) CARD(LWORD(CARD):)=' + decimal of day) -->' CALL RTNCON(CARD,LWORD(CARD)) C C End of TIME plot; begin AIRMASS plot. C CALL NEED(24) CARD(:32)='SKY brightness vs. AIRMASS, for' CARD(33:LEN+33)=BANDS(NB) CARD(LEN+34:)=DATSTR CARD(LWORD(CARD)+5:)='M = Moon up' CALL TV(CARD) CALL SPACE XLIMS(1)=1. XLIMS(2)=AIRMAX CALL PLOT(0,XLIMS,YLIMS,'L') K=0 MOONUP=0 NTWI=0 DO 30 J=JMIN,JMAX IF (KBND(J).NE.NB) GO TO 30 IF (DTYPE(J).EQ.'Y')THEN K=K+1 IF ( (DAWN.LT.DUSK .AND. 1 DJOBS(J).GT.DAWN .AND. DJOBS(J).LT.DUSK) .OR. 2 (DAWN .GT.DUSK .AND. 3 (DJOBS(J).GT.DAWN .OR. DJOBS(J).LT.DUSK) ) ) THEN C Twilight. A1='t' NTWI=NTWI+1 KX(1,K)=-1 ELSE IF ( (RISE.LT.SET .AND. 1 DJOBS(J).GT.RISE .AND. DJOBS(J).LT.SET) .OR. 2 (RISE .GT.SET .AND. 3 (DJOBS(J).GT.RISE .OR. DJOBS(J).LT.SET) ) ) THEN C Moon is up. Flag in KX. A1='M' MOONUP=MOONUP+1 C Flag for model solution. KX(1,K)=1 ELSE C Moon is down. A1='-' KX(1,K)=0 END IF SYMB(K)=A1 CALL PLOT(-1,AIRM(J),SIGNAL(J),A1) C Set KX(2,*) to NDIA, or 0 if no diaphragm: DO 28 NDIA=1,NDIAS IF (DIAPH(J).EQ.DIANAM(NDIA)) THEN KX(2,K)=NDIA NDIA2NP(NDIA,1)=NDIA2NP(NDIA,1)+1 GO TO 29 END IF 28 CONTINUE C Here only if no dias. KX(2,K)=0 C Debugging: C IF (KDIAPHRAGM.NE.-1) CALL TV('Dia.error!') 29 CONTINUE C Set KX(3,*) to J for use in later plots. KX(3,K)=J C Set X(1,*) to airmass. X(1,K)=AIRM(J) C Set X(4,*) to nearest star signal. CALL HALO(NB,DTYPE,KDIAPHRAGM,J,JBGN,JEND, X(4,K)) Y(K)=SIGNAL(J) W(K)=1./ABS(SIGNAL(J)) END IF 30 CONTINUE CALL PLOT(1,50.,-50.,' ') CALL RTNCON(' AIRMASS -->',16) C C End of AIRMASS plot; begin LUNAR ELONGATION plot. C IF(MOONUP.LT.2) GO TO 100 ELONGMAX=0. ELONGMIN=50. CALL NEED(24) CARD(:41)='SKY brightness vs. LUNAR ELONGATION, for' CARD(42:LEN+42)=BANDS(NB) CARD(LEN+43:)=DATSTR CARD(LWORD(CARD)+5:)='M = Moon up' CALL TV(CARD) CALL SPACE XLIMS(1)=0. XLIMS(2)=180. CALL PLOT(0,XLIMS,YLIMS,'L') K=0 DO 40 J=JMIN,JMAX IF (KBND(J).NE.NB) GO TO 40 IF (DTYPE(J).EQ.'Y')THEN K=K+1 T=(DTZERO(NIGHT)-51544.5D0+DJOBS(J))/36525. UTRAD=MOD(DJOBS(J),1.)*TWOPI TLST=STUTZ+(UTRAD*1.00273791)+ELONG C MOON uses T for lunar position, TLST for topocentric ephemeris. CALL MOON(T, TLST) C MOON returns lunar coords. in RASUN, DESUN in /CSUN/. HAMOON=TLST-RASUN COSDEL=COS(DESUN) SINDEL=SIN(DESUN) COSHA=COS(HAMOON) XM= -COSDEL*SIN(HAMOON) YM=COSPHI*SINDEL-SINPHI*COSDEL*COSHA ZM=SINPHI*SINDEL+COSPHI*COSDEL*COSHA IF (ZM.GT.ALTS(1)) THEN C Set X(2,*) to LUNAR airmass. (ZM is cosZ(Moon).) X(2,K)=(.0096467+ZM*(.148386+ZM*1.0024342))/ 1 (.000303978+ZM*(.0102963+ZM*(.149864+ZM))) ELSE C Moon below horizon. X(2,K)=50. C WRITE(CARD,*)'ZM=',ZM C IF(KX(1,K).EQ.1)CALL TV(CARD) IF(KX(1,K).EQ.1) SYMB(K)='m' END IF AZMOON=ATAN2(YM,XM) DELAZ=AZMOON-AZIM(J) C approximate cosz(sky) as 1/airmass. COSEL=ZM/AIRM(J) + 1 SQRT((1.-ZM*ZM)*(1.-1./AIRM(J)**2))*COS(DELAZ) C Set X(3,*) to lunar elongation. X(3,K)=ACOS(COSEL) IF (KX(1,K).EQ.1) THEN ELONGMAX=MAX(ELONGMAX,X(3,K)) ELONGMIN=MIN(ELONGMIN,X(3,K)) END IF XS(K)=ACOS(COSEL)/DEGRAD YS(K)=SIGNAL(J) END IF 40 CONTINUE CALL PLOT(0,0.,0.,'D') CALL PLOT(K,XS,YS,SYMB) CALL PLOT(0,0.,0.,'S') CALL RTNCON(' LUNAR ELONGATION (deg.) -->',32) C C End of plots. C 100 CALL SPACE CALL TVN('Do you want to interpolate in TIME alone, or') CALL TVN('do you want to try fitting sky with a MODEL?') CALL TV('(Enter T for Time; M to Model the sky; or)') CALL ASKN('( R to re-display the plots.)',A1) CALL SPACE C 102 IF (A1.EQ.'T') THEN K=0 DO 105 J=JMIN,JMAX C C restore sky data for current night and passband. IF (KBND(J).NE.NB) GO TO 105 IF (DTYPE(J).EQ.'Y')THEN C add to XS, YS. K=K+1 XS(K)=DJOBS(J) YS(K)=SIGNAL(J) END IF C 105 CONTINUE C 110 CALL TV('Enter C to replace sky data with a Constant,') CALL TVN(' L to smooth with a Linear fit,') CALL TVN(' N to use the sky Nearest the star,') CALL TVN(' P to use the sky Preceding the star,') CALL TVN(' F to use the sky Following the star,') CALL TVN(' R to Re-display plots,') CALL TVN(' or Q to Quit.') CALL ASKN('?',A1) IF (A1.EQ.'C') THEN C be lazy: fit line, but force it horizontal. CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE) SLOPE=0. PG(2)=0. NIX=1 IX(1)=2 ELSE IF (A1.EQ.'L') THEN IF (K.LT.6) THEN CALL TV('Not enough data to fit a line safely.') CALL TVN('Try a constant instead.') GO TO 110 END IF CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE) PG(2)=SLOPE NIX=0 ELSE C for neighboring-sample subtraction: C IF (A1.EQ.'N') THEN C use NEAREST neighbor. CALL TV('CAUTION: This method is not robust!') CALL SKYNBR(NB,DTYPE,KDIAPHRAGM,JBGN,JEND) ELSE IF (A1.EQ.'P') THEN C SKY comes first. CALL TV('CAUTION: This method is not robust!') CALL SKYADJ(NB,DTYPE,KDIAPHRAGM,JBGN,JEND,NERROR) IF (NERROR.GT.0) GO TO 2 ELSE IF (A1.EQ.'F') THEN C STAR comes first, so scan in reverse order. CALL TV('CAUTION: This method is not robust!') CALL SKYADJ(NB,DTYPE,KDIAPHRAGM,JEND,JBGN,NERROR) IF (NERROR.GT.0) GO TO 2 ELSE IF (A1.EQ.'R') THEN GO TO 2 ELSE CALL TV('Usually N is your best choice; but') CALL TVN('if you have been careful to take sky') CALL TVN('ALWAYS before or after every star, you') CALL TVN('can choose P for sky-then-star, or') CALL TVN(' F for star-then-sky.') GO TO 110 END IF C skip residual plots if sky is subtracted exactly. CALL TV('(Sky samples were subtracted exactly; so') CALL TVN(' there are no residuals to plot.)') GO TO 999 END IF C C here when linear model was used. C IPR=0 PG(1)=YBAR C PG(2) was set above. DO 120 J=1,K X(1,J)=XS(J)-XBAR Y(J)=YS(J) W(J)=1. 120 CONTINUE C refine fit: C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST) CALL DLSQ(YPLIN,K,2, 1, 2, NIX, 1, 0, 1.E-3) SLOPE=P(2) YBAR=P(1) C YMAX=-3.E33 YMIN=3.E33 C DO 130 J=JBGN,JEND IF (KBND(J).NE.NB) GO TO 130 C subtract fitted sky: SIGNAL(J)=SIGNAL(J)-((DJOBS(J)-XBAR)*SLOPE+YBAR) IF (DTYPE(J).EQ.'Y') THEN YMAX=MAX(YMAX,SIGNAL(J)) YMIN=MIN(YMIN,SIGNAL(J)) END IF 130 CONTINUE C C plot sky resids. C CALL NEED(24) CARD(:26)='Sky resids. vs. TIME for' CARD(27:)=BANDS(NB) CALL TV(CARD) C set plot limits.... XLIMS(1)=DJOBS(JMIN) XLIMS(2)=DJOBS(JMAX) YLIMS(1)=YMIN YLIMS(2)=YMAX CALL PLOT(0,XLIMS,YLIMS,'L') DO 140 J=JBGN,JEND IF (KBND(J).NE.NB) GO TO 140 CALL PLOT(-1,DJOBS(J),SIGNAL(J),'*') 140 CONTINUE C force plot: CALL PLOT(1,0.,0.,'+') CALL RTNCON(' Time (decimal days) -->',28) GO TO 999 C ELSE IF (A1.EQ.'M') THEN C go on. C ELSE IF (A1.EQ.'R') THEN GO TO 2 C ELSE IF (A1.EQ.'H') THEN CALL SPACE2 CALL TV('Reply M to detect bad sky measurements.') CALL TV('Reply T to use a neighboring sky sample,') CALL TVN('either before, after, or nearest the star.') CALL SPACE2 GO TO 100 C ELSE GO TO 100 END IF ELSE IF (K.EQ.1) THEN CALL TV('Only 1 datum available for this night!') CALL SPACE2 DO 150 J=JBGN,JEND C subtract sky: IF (KBND(J).EQ.NB) 1 SIGNAL(J)=SIGNAL(J)-SIGMIN 150 CONTINUE C skip plot. GO TO 999 ELSE CALL TV('NO SKY DATA AVAILABLE for this night') Call ASK('Do you want to continue?',A1) IF (A1.EQ.'Y') THEN RETURN ELSE CALL STSEPI END IF END IF C C C Now ***** MODEL ***** and interpolate sky. C IPR=0 NPTS=K C remove const. term: NIX=1 IX(1)=1 PG(1)=0. C set up diaphragm factors: NPARMS=11 NN=0 NREFDIA=0 DO 157 K=12,11+NDIAS PG(K)=1. NDIA=K-11 WRITE(CARD,*)NDIA2NP(NDIA,1),' data in dia.',DIANAM(NDIA) CALL TV(CARD) C Set NPARMS to largest one actually used. IF (NDIA2NP(NDIA,1).GT.0) NPARMS=K C Find leading diaphragm... IF (NDIA2NP(NDIA,1).GT.NN) THEN C flag NDIA as possible NREFDIA. NN=NDIA2NP(NDIA,1) NREFDIA=NDIA END IF 157 CONTINUE C C C Plot star haloes: K=0 DO 160 J=1,NPTS IF (X(4,J).GT.0.) THEN K=K+1 XS(K)=X(4,J) YS(K)=Y(J) END IF 160 CONTINUE CALL ROBLIN(XS,YS,K, XBAR, YBAR, SLOPE) C Fraction of star included: IF (SLOPE.GT.0. .AND. SLOPE.LT.0.01) THEN CALL SORT2(XS,YS,K) ZM=XS(K) C ZM is max. value of XS. CALL NEED(24) CALL CENTER('Sky vs. corresponding Star signal:') CALL TVN(' Sky') CALL PLOT(-K,XS,YS,' ') CALL PLOT(0,0.,0.,'O') DO 165 J=1,20 XM=J*ZM/20. YM=YBAR+SLOPE*(XM-XBAR) CALL PLOT(-1,XM,YM,'+') 165 CONTINUE CALL PLOT(0,0.,1.,'O') DO 166 J=1,NPTS IF (X(4,J).EQ.0.) GO TO 166 CALL PLOT(-1,X(4,J),Y(J),SYMB(J)) 166 CONTINUE CALL PLOT(1,0.,0.,'+') CARD=' Star signal --> (Sky is 0.xxxx of Star)' WRITE (CARD(32:38),'(F7.4)') SLOPE CALL RTNCON(CARD,49) ELSE CALL TV('Contribution of star to sky measures not detected.') SLOPE=0. END IF C C Set default order: DO 167 J=1,10 LORDER(J)=J 167 CONTINUE C C estimate PG's: C PG(11)=SLOPE C C Start for dark sky: C Fit rough line to dark sky: K=0 DO 170 J=1,NPTS IF (KX(1,J).EQ.0) THEN C Add to new list: K=K+1 XS(K)=X(1,J) YS(K)=Y(J) - PG(11)*X(4,J) C set initial wt.: W(J)=1. ELSE C zero weights if Moon is up, or twilight: W(J)=0. END IF 170 CONTINUE NN=K C IF (NPTS-MOONUP-NTWI.GT.6 .AND. AIRMAX.GT.1.5) THEN C Debugging: C CALL TV('DATA:') C CALL PLOT(-NN,XS,YS,'*') C NN data now set for moonless part. C OLDVAR=3.E33 NFIT=1 C Set order of relaxation: LORDER(3)=4 LORDER(4)=5 LORDER(5)=3 C C Fit line to data vs. airmass: CALL ROBLIN(XS,YS,NN, XBAR,YBAR,SLOPE) C Debugging: C DO 174 K=1,NN C CALL PLOT(-1, XS(K), YBAR+SLOPE*(XS(K)-XBAR),'-') C 174 CONTINUE C CALL PLOT(1,0.,0.,'+') C CALL RTNCON(' AIRMASS -->',13) PG(2)=YBAR/XBAR C IF (YBAR.LE.XBAR*SLOPE) THEN C Negative or zero intercept; form can't be fit. C CALL TV('Negative intercept.') PG(3)=0.D0 PG(4)=0.D0 PG(5)=0.D0 ELSE C Positive intercept; OK. IF (SLOPE.GT.0.5*YBAR/XBAR) THEN C P(4) is small. C CALL TV('P(4) small.') PG(4)=(SQRT(YBAR/(XBAR*SLOPE))-1.)/XBAR PG(5)=PG(4)**2 ELSE C P(4) is big. C CALL TV('P(4) big; using P(5).') PG(4)=1./XBAR PG(5)=(YBAR/XBAR-SLOPE)/((YBAR+SLOPE*XBAR)*XBAR) PG(2)=(YBAR/XBAR+YBAR*PG(4)+YBAR*XBAR*PG(5)) LORDER(3)=5 LORDER(4)=4 END IF END IF C Override trial line with a priori values: PG(4)=EXTIN(NB) PG(5)=0.5*PG(4)*PG(4) PG(2)=YBAR/(1.+(PG(4)+PG(5)*XBAR)*XBAR) ELSE C Not enough data or range to fit. CALL TV('Can''t fit dark sky.') PG(2)=YMIN*2. PG(4)=EXTIN(NB) PG(5)=0.5*PG(4)*PG(4) END IF PG(3)=PG(5)*PG(2)*0.1 ZENITH=PG(1)+(PG(2)+PG(3))/(1.+PG(4)+PG(5)) IF (NPTS-MOONUP-NTWI.LT.6) GO TO 200 C C PG's now set. C C hold all fixed but P(2): DO 180 J=2,11 CALL FIXP(J,PG(J),NIX) 180 CONTINUE CALL UNFIXP(2,NIX) C C Relax successive parameters: DO 190 L=3,6 CARD='Solving for x sky parameters...' C IF(L.EQ.3) CARD(24:24)=' ' WRITE(CARD(12:13),'(I2)') L-1 CALL TV(CARD) IF (NPTS+NIX.LT.NPARMS+4) GO TO 191 C C Fix all sky parms.: DO 181 K=1,11 IX(K)=K 181 CONTINUE NIX=11 IF (NDIAS.GT.0) THEN C Fix reference diaphragm: K=NREFDIA+11 CALL FIXP(K,PG(K),NIX) C Debugging: C WRITE(CARD,*)'Holding PG(',K,') at',PG(K) C CALL TV(CARD) END IF C C See if we can release P(11), the star fraction: IF (NPTS.GT.2*MOONUP) CALL UNFIXP(11,NIX) C C First, vary ONLY the new parameter. CALL UNFIXP(LORDER(L-1),NIX) C C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST) CALL DLSQ(YPSKY,NPTS,NPARMS, 1, NPARMS,NIX, 4, 2, 1.E-3) C C Next, see if result is acceptable: IF (WVAR.GT.0.D0) PG(LORDER(L-1))=P(LORDER(L-1)) C Check each case individually: IF (P(2).LT.0.D0) CALL FIXP(2,0.,NIX) IF (P(3).LT.0.D0) CALL FIXP(3,0.,NIX) IF (P(5).LT.0.D0) CALL FIXP(5,0.,NIX) IF (P(11).LT.0.D0) CALL FIXP(11,0.,NIX) C DO 182 K=3,L CALL UNFIXP(LORDER(K-1),NIX) 182 CONTINUE C C Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST) CALL DLSQ(YPSKY,NPTS,NPARMS, 1, NPARMS,NIX, 4, 2, 1.E-3) C C Give up if solution failed: IF (WVAR.EQ.0.D0 .OR. WVAR.GT.1.2*OLDVAR) GO TO 189 C C Debugging: C CALL TV('DATA:') C CALL PLOT(-NN,XS,YS,'*') C DO 184 K=1,NPTS C remember that YC is DP... C IF (W(K).GT.0.) CALL PLOT(-1,X(1,K),REAL(YC(K)),'C') C 184 CONTINUE C CALL PLOT(1,0.,0.,'+') C CALL RTNCON(' AIRMASS -->',13) C C Check solution for plausibility: IF (P(2).LE.0.D0 .OR. P(3).LT.0.D0 .OR. P(5).LT.0.D0 .OR. 1 (P(4).LT.0.D0 .AND. P(4)*P(4).GE.4.D0*P(5)) .OR. 2 (P(4).GT.100.D0 .OR. P(5).GE.100.D0) .OR. 3 (P(11).GT.0.01D0 .OR. P(11).LT.0.D0) ) THEN C Trouble. GO TO 189 ELSE C OK. END IF C C Adjust weights of moonless data: ZENITH=PG(1)+(PG(2)+PG(3))/(1.+PG(4)+PG(5)) DO 185 K=1,NPTS IF (KX(1,K).EQ.0) W(K)=EXPTIM(KX(3,K))/(ABS(YC(K)) + ZENITH) 185 CONTINUE C Update parameters: DO 187 K=1,5 PG(K)=P(K) ERRSAVE(K)=SP(K) C WRITE(CARD,'(''P('',I2,'') ='',G9.2)') K,P(K) C CALL TV(CARD) C WRITE(CARD,'(4X,''+/-'',G9.2)') SP(K) C CALL TVN(CARD) 187 CONTINUE C Star wings: PG(11)=P(11) C IF (SP(11).GT.0.) THEN C WRITE(CARD,'(''P(11) ='',G9.2)') P(11) C CALL TV(CARD) C WRITE(CARD,'(4X,''+/-'',G9.2)') SP(11) C CALL TVN(CARD) C END IF C Diaphragms: DO 188 K=12,NPARMS IF (P(K).GT.0.) PG(K)=P(K) 188 CONTINUE NFIT=NPARMS-NIX WRITE(CARD,'(I5,A)') NFIT,'-parameter fit accepted.' CALL TV(CARD) OLDVAR=WVAR C Release next parameter and try again. 189 CALL UNFIXP(LORDER(L),NIX) 190 CONTINUE C 191 CONTINUE C C Plot fit: CALL SPACE CALL NEED(24) WRITE (CARD,'(A,I3,A,I4,2A)') 1 'Fit',NFIT,' parameters to',NN,' dark-sky data in ',BANDS(NB) CALL CENTER(CARD) CALL TVN(' Sky') ZENITH=PG(1)+(PG(2)+PG(3))/(1.+PG(4)+PG(5)) C Reserve space: XLIMS(1)=1.0 XLIMS(2)=AIRMAX YLIMS(1)=0.7*ZENITH XM=SQRT(AIRMAX) YLIMS(2)=1.4*MAX(ZENITH, PG(1)+((PG(2)+AIRMAX*PG(3))*AIRMAX)/ 1 (1.+AIRMAX*(PG(4)+AIRMAX*PG(5))), 2 PG(1)+((PG(2)+XM*PG(3))*XM)/(1.+XM*(PG(4)+XM*PG(5)))) CALL PLOT(0,XLIMS,YLIMS,'L') CALL PLOT(0,0.,0.,'O') C Show fitted fcn.: DO 194 K=1,40 XS(1)=1.+(AIRMAX-1.)*.025*K YS(1)=PG(1)+((PG(2)+XS(1)*PG(3))*XS(1))/ 1 (1.+XS(1)*(PG(4)+XS(1)*PG(5))) CALL PLOT(-1,XS,YS,'+') 194 CONTINUE CALL PLOT(0,0.,1.,'O') C Show data: DO 195 K=1,NPTS IF(KX(1,K).EQ.0) CALL PLOT(-1,X(1,K),Y(K)-PG(11)*X(4,K), '*') 195 CONTINUE CALL PLOT(1,0.,0.,'+') CALL RTNCON(' Airmass -->',14) C C Start for MOONLIGHT: C C Rough estimates for large-airmass parameters: 200 PG(9)=0.35*EXTIN(NB) PG(10)=EXP(-EXTIN(NB)) C C Fit rough fcn. to lunar aureole: K=0 DO 210 J=1,NPTS IF (KX(1,J).EQ.1) THEN C Add to new list: K=K+1 XM=X(1,J)+X(2,J) ZM=X(1,J)*X(2,J) XS(K)=X(3,J) C remove dark-sky component... YS(K)=(Y(J) - (PG(1)+((PG(2)+X(1,J)*PG(3))*X(1,J))/ 1 (1.+X(1,J)*(PG(4)+X(1,J)*PG(5)))) C ... remove star halo ... 2 - PG(11)*X(4,J) C ...and scale to unit airmass. 3 ) * X(3,J)/(X(1,J)*(EXP(-PG(9)*XM)+PG(10)/ZM)) C partial allowance for photon noise: W(J)=EXPTIM(KX(3,J))/(ABS(Y(J))+ZENITH) ELSE C Keep old wt. for dark sky. END IF 210 CONTINUE C IF (MOONUP.GT.6 .AND. ELONGMAX-ELONGMIN.GT.0.2) THEN C Sort on elongation & use bottom of list for aureole: CALL SORT2(XS,YS,K) DO 211 J=1,K C Select points in aureole: IF (XS(J).GT.1.) GO TO 212 211 CONTINUE C Whole set is in aureole. Use 1st half: J=NPTS/2 212 IF (J.EQ.1) THEN C No aureole data. Pick something: J=3 END IF C Debugging: C WRITE(CARD,*) J,' points in initial fit.' C CALL TV(card) C First J now suitable to estimate PG(6). CALL ROBLIN(XS,YS,J, XBAR,YBAR,SLOPE) PG(6)=YBAR-XBAR*SLOPE PG(7)=SLOPE IF (PG(6).LT.0.D0) THEN C Aureole not detectable. PG(6)=0.D0 PG(7)=YBAR/XBAR END IF C Now move down *last* third... DO 216 J=1,K/3 XS(J)=XS(K-K/3+J) C ...and divide by elongation. YS(J)=YS(K-K/3+J)/XS(J) 216 CONTINUE CALL ROBLIN(XS,YS,K/3, XBAR,YBAR,SLOPE) PG(8)=MAX(SLOPE,0.) ELSE C Not enough data to fit with Moon. IF (K.GT.0) THEN C Use median Y: CALL SORT2(YS,XS,K) PG(7)=(YS((K+1)/2)+YS(K/2+1))*0.5 ELSE PG(7)=0. END IF PG(6)=0. PG(8)=0. END IF C C PG's now set for moonlight terms. C IF (ELONGMAX-ELONGMIN.LT.2.) THEN C Rearrange order of solving: IF (ELONGMIN.GT.1.5) THEN LORDER(6)=8 LORDER(8)=10 LORDER(10)=6 IF (ELONGMIN.GT.2.) THEN LORDER(7)=9 LORDER(9)=7 END IF ELSE IF (ELONGMIN.GT.0.5) THEN LORDER(6)=7 LORDER(7)=6 IF (ELONGMAX.LT.1.) THEN LORDER(7)=9 LORDER(8)=10 LORDER(9)=6 LORDER(10)=8 END IF ELSE IF (ELONGMAX.LT.1.0) THEN LORDER(7)=9 LORDER(8)=7 LORDER(9)=10 LORDER(10)=8 END IF END IF C OLDVAR=3.E33 NFIT=0 IF (MOONUP.LT.3) GO TO 241 DO 240 J=6,10 C Relax parameters one by one: DO 221 K=1,11 IX(K)=K 221 CONTINUE NIX=11 IF (NDIAS.GT.0) THEN C Fix reference diaphragm: K=NREFDIA+11 CALL FIXP(K,PG(K),NIX) END IF C C Release P(11), the star fraction: CALL UNFIXP(11,NIX) C CARD='Solving for x moonlight parameters...' IF(J.EQ.6) CARD(24:24)=' ' WRITE(CARD(12:13),'(I2)') J-5 CALL TV(CARD) C C First, vary ONLY the new parameter. CALL UNFIXP(LORDER(J),NIX) IF (MOONUP.LT.NPARMS+4-NIX) GO TO 241 C C subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST) CALL DLSQ(YPSKY,NPTS,NPARMS, 1,NPARMS,NIX, 4, 2, 1.E-3) C C Next, see if result is acceptable: IF (WVAR.GT.0.D0) PG(LORDER(J))=P(LORDER(J)) C Check each case individually: IF (P(6).LT.0.D0) CALL FIXP(6,0.,NIX) IF (P(8).LT.0.D0) CALL FIXP(8,0.,NIX) IF (P(9).LT.0.D0) CALL FIXP(9,0.,NIX) IF (P(10).LT.0.D0) CALL FIXP(10,0.,NIX) IF (P(9).GT.9.D0 .AND. P(10).EQ.0.D0) CALL FIXP(9,9.,NIX) IF (P(9).GT.10.D0) CALL FIXP(9,1.,NIX) IF (P(11).LT.0.D0) CALL FIXP(11,0.,NIX) C C Then, vary parameters 6 to J. DO 222 K=6,J CALL UNFIXP(LORDER(K),NIX) 222 CONTINUE IF (MOONUP.LT.NPARMS+4-NIX) GO TO 241 C C subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST) CALL DLSQ(YPSKY,NPTS,NPARMS, 1,NPARMS,NIX, 4, 2, 1.E-3) C C Re-extract data for plots: YMAX=-3.E33 YMIN=3.E33 K=0 DO 226 L=1,NPTS IF (KX(1,L).EQ.1) THEN C Add to new list: K=K+1 XS(K)=X(3,L) XM=X(1,L)+X(2,L) ZM=X(1,L)*X(2,L) C remove dark-sky component... YS(K)=(Y(L) - (P(1)+((P(2)+X(1,L)*P(3))*X(1,L))/ 1 (1.+X(1,L)*(P(4)+X(1,L)*P(5)))) C ...subtract star halo, and scale to unit airmass. 2 - PG(11)*X(4,L) ) / (X(1,L)*(EXP(-P(9)*XM) + P(10)/ZM)) IF(NDIAS.GT.1) YS(K)=YS(K)/P(KX(2,L)+11) YMAX=MAX(YMAX,YS(K)) YMIN=MIN(YMIN,YS(K)) END IF 226 CONTINUE NN=K C C Skip plotting: C IF (NIX.GT.0) GO TO 231 C C Fix up plot limits: C XLIMS(1)=0.1 XLIMS(2)=ELONGMAX C Adjust for possible outliers: CALL SORT2(YS,XS,K) C (Note reversal of Y and X.) YM=(YS((K+1)/2)+YS(K/2+1))*0.5 C YM is median. L=(K+2)/4 ZM=YS(K-L)-YS(L) C ZM is interquartile range. YMAX=MIN(YMAX,YM+4.*ZM) YMIN=MAX(YMIN,YM-4.*ZM) YLIMS(1)=YMIN YLIMS(2)=YMAX C C Show result: CALL SPACE CALL NEED(24) WRITE(CARD,'(I5,A,I4,1X,2A)') 1 NPARMS-NIX,'-parameter fit ( + + + ) to',MOONUP, 2 BANDS(NB)(:LWORD(BANDS(NB))),' data ( * ) with Moon' CALL CENTER(CARD) CALL TVN('Normalized Sky') CALL PLOT(0,0.,0.,'O') C Reserve space: CALL PLOT(0,XLIMS,YLIMS,'L') IF (YLIMS(1).LT.0.) CALL XAXIS(XLIMS) C Plot fit to aureole terms: DO 230 K=1,40 XM=0.1+(.02+.003*K)*K YM=(P(6)/XM + P(7) + P(8)*XM) CALL PLOT(-1,XM,YM,'+') 230 CONTINUE CALL PLOT(0,0.,1.,'O') CALL PLOT(NN,XS,YS,'*') CALL RTNCON(' LUNAR ELONGATION (rad.) -->',32) C C Give up if solution failed: 231 IF (WVAR.EQ.0.D0 .OR. WVAR.GT.1.2*OLDVAR) GO TO 240 C Tests for plausibility: IF (P(6).LT.0.D0 .OR. P(8).LT.0.D0 .OR. 1 P(9).LT.0.D0 .OR. P(10).LT.0.D0 .OR. 2 (P(9).GT.10.D0 .AND. P(10).EQ.0.D0) .OR. 3 (P(11).GT.0.01D0 .OR. P(11).LT.0.D0) ) GO TO 240 C C Result accepted. C C Revise weights: DO 235 K=1,NPTS C partial allowance for photon noise: IF(KX(1,K).EQ.1) W(K)=EXPTIM(KX(3,K))/(ABS(Y(K))+ZENITH) 235 CONTINUE C Update parameters: DO 237 K=6,10 PG(K)=P(K) ERRSAVE(K)=SP(K) C WRITE(CARD,'(''P('',I2,'') ='',G9.2)') K,P(K) C CALL TV(CARD) C WRITE(CARD,'(4X,''+/-'',G9.2)') SP(K) C CALL TVN(CARD) 237 CONTINUE C Wings of star image: IF (P(11).GT.0.D0) PG(11)=P(11) IF (SP(11).GT.0.) THEN C WRITE(CARD,'(''P(11) ='',G9.2)') P(11) C CALL TV(CARD) C WRITE(CARD,'(4X,''+/-'',G9.2)') SP(11) C CALL TVN(CARD) END IF C Diaphragms: DO 238 K=12,NPARMS IF (P(K).GT.0.) PG(K)=P(K) 238 CONTINUE C NFIT=NPARMS-NIX WRITE(CARD,'(I5,A)') NFIT,'-parameter fit accepted.' CALL TV(CARD) OLDVAR=WVAR C 240 CONTINUE C C Re-extract data for final plot: 241 YMAX=-3.E33 YMIN=3.E33 K=0 DO 244 L=1,NPTS IF (KX(1,L).EQ.1) THEN C Add to new list: K=K+1 XS(K)=X(3,L) XM=X(1,L)+X(2,L) ZM=X(1,L)*X(2,L) C remove dark-sky component... YS(K)=(Y(L) - (PG(1)+((PG(2)+X(1,L)*PG(3))*X(1,L))/ 1 (1.+X(1,L)*(PG(4)+X(1,L)*PG(5)))) C ...subtract star halo, and scale to unit airmass. 2 - PG(11)*X(4,L) ) / (X(1,L)*(EXP(-PG(9)*XM) + PG(10)/ZM)) YMAX=MAX(YMAX,YS(K)) YMIN=MIN(YMIN,YS(K)) END IF 244 CONTINUE C C Fix up plot limits: XLIMS(1)=ELONGMIN XLIMS(2)=ELONGMAX IF (K.EQ.0) GO TO 300 C Adjust for possible outliers: CALL SORT2(YS,XS,K) C (Note reversal of Y and X.) YM=(YS((K+1)/2)+YS(K/2+1))*0.5 C YM is median. L=(K+2)/4 ZM=YS(K-L)-YS(L) C ZM is interquartile range. YMAX=MIN(YMAX,YM+3.*ZM) YMIN=MAX(YMIN,YM-3.*ZM) IF (YMAX.EQ.YMIN) GO TO 300 YLIMS(1)=YMIN YLIMS(2)=YMAX C C Show result: CALL SPACE CALL NEED(24) WRITE(CARD,'(I5,A,I4,1X,2A)') 1 NFIT,'-parameter fit ( + + + ) to',MOONUP, 2 BANDS(NB)(:LWORD(BANDS(NB))),' data ( * ) with Moon' CALL CENTER(CARD) CALL TVN('Normalized Sky') CALL PLOT(0,0.,0.,'O') C Reserve space: CALL PLOT(0,XLIMS,YLIMS,'L') IF (YLIMS(1).LT.0.) CALL XAXIS(XLIMS) C Plot fit to aureole terms: DO 250 K=1,40 XM=ELONGMIN+(ELONGMAX-ELONGMIN)*(.02+.003*K)*K YM=(PG(6)/XM + PG(7) + PG(8)*XM) CALL PLOT(-1,XM,YM,'+') 250 CONTINUE CALL PLOT(0,0.,1.,'O') CALL PLOT(NN,XS,YS,'*') CALL RTNCON(' LUNAR ELONGATION (rad.) -->',32) C 300 CONTINUE C C C Now display fitted equation: C DO 310 J=1,10 SP(J)=ERRSAVE(J) 310 CONTINUE CALL SKYEQN(PG,SP,CARD) C C C Don't forget diaphragm scaling: C IF (NPARMS.GT.11) THEN CALL SPACE CALL TV('Diaphragm scaling:') END IF DO 320 J=12,NPARMS WRITE(CARD,'(A,F7.3,'' +/-'',F7.3)') DIANAM(J-11),P(J),SP(J) CALL TV(CARD) 320 CONTINUE C C C Now display overall fit: C C Hold all fixed; evaluate YC's: DO 355 K=1,NPARMS IX(K)=K 355 CONTINUE NIX=NPARMS C C subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST) CALL DLSQ(YPSKY,NPTS,NPARMS, 1, NPARMS,NIX, 4, 2, 1.E-1) C CALL SPACE CALL NEED(24) CARD='Measured vs. model ( + + ) sky brightness in '//BANDS(NB) CARD(LWORD(CARD)+1:)=':' CALL CENTER(CARD) CALL TVN('Observed sky') XLIMS(1)=3.E33 YLIMS(1)=0. YLIMS(2)=0. K=0 DO 360 J=1,NPTS C omit twilight. IF (KX(1,J).LT.0) GO TO 360 K=K+1 XS(K)=YC(J) YS(K)=Y(J) XLIMS(1)=MIN(XLIMS(1),XS(K)) IF (KX(1,J).EQ.0) THEN SYMB(K)='*' ELSE IF (KX(1,J).EQ.1) THEN SYMB(K)='M' C ELSE IF (KX(1,J).EQ.-1) THEN C SYMB(K)='t' ELSE SYMB(K)='?' END IF 360 CONTINUE CALL SORT2(XS,YS,K) IF (XS(K).LT.2.*XS(K-1)) THEN XLIMS(2)=XS(K) ELSE XLIMS(2)=XS(K-1) END IF YLIMS(1)=0.8*XLIMS(1) YLIMS(2)=1.2*XLIMS(2) CALL PLOT(0,XLIMS,YLIMS,'L') C Draw unit diagonal: DO 370 J=1,20 XM=J*XLIMS(2)/20. CALL PLOT(-1,XM,XM,'+') 370 CONTINUE CALL PLOT(0,0.,0.,'D') CALL PLOT(K,XS,YS,SYMB) CALL RTNCON(' Calculated sky -->',23) C C Plot star halo resids.: C CALL SPACE CALL NEED(24) CARD='Residuals vs. Star signal in '//BANDS(NB) CARD(LWORD(CARD)+1:)=':' CALL CENTER(CARD) CALL TVN(' Residual') XLIMS(1)=0. YLIMS(1)=0. YLIMS(1)=0. YLIMS(2)=0. K=0 DO 380 J=1,NPTS IF (KX(1,J).LT.0) GO TO 380 K=K+1 XS(K)=X(4,J) YS(K)=Y(J)-YC(J) XLIMS(2)=MAX(XLIMS(2),XS(K)) YLIMS(1)=MIN(YLIMS(1),YS(K)) YLIMS(2)=MAX(YLIMS(2),YS(K)) 380 CONTINUE CALL PLOT(0,XLIMS,YLIMS,'L') CALL XAXIS(XLIMS) CALL PLOT(K,XS,YS,SYMB) CALL PLOT(0,0.,0.,'S') CARD=' Star signal --> (Sky is 0.xxxx of Star)' WRITE (CARD(32:38),'(F7.4)') PG(11) CALL RTNCON(CARD,49) C C Plot resids. vs. time: C CALL SPACE CALL NEED(24) CARD='O/C Ratio vs. Time in '//BANDS(NB) CARD(LWORD(CARD)+1:)=':' CALL CENTER(CARD) CALL TVN(' Ratio') XLIMS(1)=DJOBS(JDUSK) XLIMS(2)=DJOBS(JDAWN) YLIMS(1)=0.6 YLIMS(2)=1.5 CALL PLOT(0,XLIMS,YLIMS,'L') C Draw "axis" at 1.0: DO 388 J=1,80 XM=XLIMS(1)+(XLIMS(2)-XLIMS(1))*J/80. CALL PLOT(-1,XM,1.,'-') 388 CONTINUE DO 390 J=1,NPTS IF (KX(1,J).LT.0) GO TO 390 IF (YC(J).EQ.0.D0) YC(J)=1.E-6 CALL PLOT(-1,DJOBS(KX(3,J)),Y(J)/REAL(YC(J)),SYMB(J)) 390 CONTINUE C force plot: CALL PLOT(1,10.,0.,'+') WRITE(CARD,'(A,F7.0)')' TIME (MJD =',DTZERO(NIGHT) CARD(LWORD(CARD):)=' + decimal of day) -->' CALL RTNCON(CARD,LWORD(CARD)) C C C List outliers: C CALL SPACE2 CALL NEED(10) CARD='Sky anomalies in '//BANDS(NB) WRITE(CARD(LWORD(CARD)+2:),'(''on MJD'',F7.0)') 1 DTZERO(NIGHT)+INT(DJOBS((JBGN+JEND)/2)) CALL TV(CARD) CARD='DECIMAL SKY LUNAR LUNAR STELLAR' CARD(52:)='SKY SKY O/C' CALL TV(CARD) CARD=' DAY AIRMASS AIRMASS ELONGATION SIGNAL OBSERVED' CARD(59:)='CALCULATED RATIO' CALL TVN(CARD) WRITE(CARD,'(''------'',7(3X,''-------''))') CALL TVN(CARD) NN=0 DO 450 J=1,NPTS IF (KX(1,J).GE.0 .AND. 1 W(J)*(Y(J)-REAL(YC(J)))**2 .GT. SCALE*36.) THEN WRITE(CARD,'(F6.4,F9.3,2F10.3,1PE12.2,0P2F10.3,F8.3)') 1 DJOBS(KX(3,J)),(X(K,J),K=1,4),Y(J),YC(J),Y(J)/YC(J) IF (KX(1,J).NE.1) THEN CARD(20:35)=' - -' END IF CALL TVN(CARD) NN=NN+1 END IF 450 CONTINUE IF (NN.EQ.0) THEN CALL TV(' (none)') ELSE WRITE(CARD,'(I3,'' anomalous sky data'')') NN IF (NN.EQ.1) CARD(22:23)='um' CALL TV(CARD) END IF C C C Ask for approval: C 500 CARD= 1 'Do you want to subtract the fitted model from stellar data?' CALL ASK(CARD,A1) IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN C Go on. ELSE IF (A1.EQ.'N') THEN C Take short-cut to 110. SHORT=.TRUE. GO TO 3 ELSE IF (A1.EQ.'R') THEN GO TO 2 ELSE IF (A1.EQ.'H') THEN IF (NPTS.GT.20*NN .AND. SCALE.LT.0.1*ZENITH) THEN C Clearly a good fit. CALL TV(' YES, that''s a good fit.') ELSE IF (NPTS.LT.10*NN .OR. SCALE.GT.0.2*ZENITH) THEN C Clearly a bad fit. CALL TV(' NO, that''s a POOR fit.') ELSE C Not clear. CALL SPACE WRITE(CARD,'(I5,''% rejects; typical residual is'',F5.2)') 1 (NN*100+NPTS/2)/NPTS,SCALE/ZENITH CARD(LWORD(CARD)+2:)='of zenith brightness.' CALL TV(CARD) CALL TV('Hard to say. You may want to try it both ways,') CALL TVN('and compare the results.') END IF GO TO 500 ELSE CALL SPACE2 CALL TV('Enter either YES or NO, or') CALL TVN('H for Help, or R to Re-display plots.') GO TO 500 END IF C C C Now subtract sky fit from data: C C JBGN, JEND index the first and last data of ANY kind for the night. C JMIN, JMAX index the first and last SKY data for the night & band. C C IF (NTWI.GT.0) THEN CALL TV('Using nearest-neighbor subtraction in twilight...') CALL SKYNBR(NB,DTYPE,KDIAPHRAGM,JBGN,JDUSK) CALL SKYNBR(NB,DTYPE,KDIAPHRAGM,JDAWN,JEND) END IF C C Subtract fitted fcn. from rest: CALL TV('Using fitted model...') DO 600 J=JDUSK,JDAWN IF (NB.NE.KBND(J)) GO TO 600 Z(1)=AIRM(J) IF ( (RISE.LT.SET .AND. 1 DJOBS(J).GT.RISE .AND. DJOBS(J).LT.SET) .OR. 2 (RISE .GT.SET .AND. 3 (DJOBS(J).GT.RISE .OR. DJOBS(J).LT.SET) ) ) THEN C Moon is up. KZ(1)=1 T=(DTZERO(NIGHT)-51544.5D0+DJOBS(J))/36525. UTRAD=MOD(DJOBS(J),1.)*TWOPI TLST=STUTZ+(UTRAD*1.00273791)+ELONG C MOON uses T for lunar position, TLST for topocentric ephemeris. CALL MOON(T, TLST) C MOON returns lunar coords. in RASUN, DESUN in /CSUN/. HAMOON=TLST-RASUN COSDEL=COS(DESUN) SINDEL=SIN(DESUN) COSHA=COS(HAMOON) XM= -COSDEL*SIN(HAMOON) YM=COSPHI*SINDEL-SINPHI*COSDEL*COSHA ZM=SINPHI*SINDEL+COSPHI*COSDEL*COSHA IF (ZM.GT.0.) THEN C Set Z(2) to LUNAR airmass. Z(2)=1./ZM ELSE C Moon below horizon. Z(2)=100. END IF AZMOON=ATAN2(YM,XM) DELAZ=AZMOON-AZIM(J) C approximate cosz(sky) as 1/airmass. COSEL=ZM/AIRM(J) + 1 SQRT((1.-ZM*ZM)*(1.-1./AIRM(J)**2))*COS(DELAZ) C Set Z(3) to lunar elongation. Z(3)=ACOS(COSEL) ELSE C Moon is down. KZ(1)=0 Z(2)=0. END IF C IF (DTYPE(J).EQ.'S') THEN C Star. Z(4)=SIGNAL(J) ELSE IF (DTYPE(J).EQ.'Y') THEN C Sky. CALL HALO(NB,DTYPE,KDIAPHRAGM,J,JBGN,JEND, Z(4)) ELSE C Dark, etc. GO TO 600 END IF C C Put NDIA in KZ(2): DO 580 NDIA=1,NDIAS IF (DIAPH(J).EQ.DIANAM(NDIA)) THEN KZ(2)=NDIA GO TO 581 END IF 580 CONTINUE KZ(2)=0 581 CONTINUE C C Z's and KZ(1) are now set for YPSKY: CALL YPSKY C Sky brightness now returned in YT. SIGNAL(J)=SIGNAL(J)-YT 600 CONTINUE C C C Prepare for next night: C 999 JBGN=JEND+1 C 1000 CONTINUE C END of loops over bands and nights. C RETURN END SUBROUTINE YPSKY C C Fits Y = [Y1 + Y2 + Y3 + P(1)] * P(diaphragm), where C C C 2 C P(2)*Z(1) + P(3)*Z(1) C Y1 = --------------------------- is the sky background, and C 2 C 1 + P(4)*Z(1) + P(5)*Z(1) C C C C Y2 = Z(1) * [P(6)/Z(3) + P(7) + P(8)*Z(3)] * {exp[-P(9)*v1] + P(10)/v3} C C is the scattered moonlight. C C C Y3 = Z(4) * P(11) is the part of the star included in sky measures. C C C C P(1) is a possible zero-offset, if dark was not subtracted. C C C Independent variables: C --------------------- C C Z(1) = observed airmass v1 = Z(1) + Z(2) = sum of airmasses C C Z(2) = lunar airmass v3 = Z(1)*Z(2) = product of airmasses C C Z(3) = lunar elongation C C Z(4) = nearest star signal, or zero if none found. C C C KZ(1) = Moon flag (1 if Moon up) C C KZ(2) = diaphragm number (KDIA) C C C C IMPLICIT NONE C C Includes for dlsq: INCLUDE 'MID_REL_INCL:kingfit.inc' INCLUDE 'MID_REL_INCL:dlsq.inc' C DOUBLE PRECISION DEN, V1, V3, FIRST, SECOND, FZ1, EXPFAC, Y2 C INTEGER I, KDIA C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C IF (KZ(1).LT.0) THEN C Twilight: Not in model. KIP=0 YT=0. RETURN END IF C C C 2 C P(2)*Z(1) + P(3)*Z(1) C Y1 = --------------------------- is the sky background... C 2 C 1 + P(4)*Z(1) + P(5)*Z(1) C C KIP=6 IP(1)=1 IP(2)=2 IP(3)=3 IP(4)=4 IP(5)=5 IP(6)=11 C DEN = 1.D0 + (P(4) + P(5)*Z(1))*Z(1) YT = (P(2) + P(3)*Z(1))*Z(1)/DEN C PART(1)=1.D0 PART(2)=Z(1)/DEN PART(3)=Z(1)*PART(2) PART(4)= -YT*PART(2) PART(5)= -YT*PART(3) PART(11)=Z(4) C YT = YT + P(1) + P(11)*Z(4) C IF (KZ(1).EQ.0) THEN C Don't need moonlight, so finish up: IF (KZ(2).GT.0) THEN C Scale old partials by diaphragm factor... KDIA=KZ(2)+11 DO 30 I=1,KIP PART(IP(I))=PART(IP(I))*P(KDIA) 30 CONTINUE C Scale by diaphragm factor, and add to KIP: KIP=7 IP(7)=KDIA PART(KDIA)=YT YT=P(KDIA)*YT END IF RETURN END IF C C C Here if we need to add moonlight: C C C Y2 = Z(1) * [P(6)/Z(3) + P(7) + P(8)*Z(3)] * {exp[-P(9)*v1] + P(10)/v3} C C is the scattered moonlight. C KIP=11 IP(7)=6 IP(8)=7 IP(9)=8 IP(10)=9 IP(11)=10 C V1=Z(1)+Z(2) V3=Z(1)*Z(2) C FIRST=P(6)/Z(3) + P(7) + P(8)*Z(3) FZ1=FIRST*Z(1) IF (-P(9)*V1.GT.30.D0) THEN C try to prevent floating-point exceptions. EXPFAC=EXP(30.D0) ELSE EXPFAC=EXP(-P(9)*V1) END IF SECOND=EXPFAC + P(10)/V3 C Y2 = FZ1*SECOND C PART(7)=SECOND*Z(1) PART(6)=PART(7)/Z(3) PART(8)=Z(3)*PART(7) PART(9)= -V1*EXPFAC*FZ1 PART(10)=FZ1/V3 C YT = YT + Y2 C C Diaphragm scaling: IF (KZ(2).GT.0) THEN C Scale old partials by diaphragm factor... KDIA=KZ(2)+11 DO 50 I=1,KIP PART(IP(I))=PART(IP(I))*P(KDIA) 50 CONTINUE C Scale by diaphragm factor, and add to KIP: KIP=12 IP(12)=KDIA PART(KDIA)=YT YT=P(KDIA)*YT END IF C RETURN END SUBROUTINE SKYADJ(NB,DTYPE,KDIAPHRAGM,JBGN,JEND,NERROR) C C Picks an ADJACENT sky sample, and subtracts it from star. C C IMPLICIT NONE C REAL SKY, PHNOISE, ERREST INTEGER NB, KDIAPHRAGM, JBGN, JEND, NERROR, J3, J, 1 LWORD, JJ C INCLUDE 'MID_REL_INCL:obs.inc' C INCLUDE 'MID_REL_INCL:mstars.inc' C PARAMETER (MSTARS=1650) C MSTARS is max.number of star-catalog entries. C commons for star catalog: CHARACTER *32 STARS COMMON /SCATA/ STARS(MSTARS) C CHARACTER DTYPE(MXOBS) C CHARACTER A, CARD*48, PRECEDED*8, FOLLOWING*9, DSKY*4 C C SKY=-3.E33 DSKY=' ' NERROR=0 C C Decide which way to scan: C IF (JBGN.LT.JEND) THEN J3=1 PRECEDED='preceded' FOLLOWING='FOLLOWING' ELSE J3=-1 PRECEDED='followed' FOLLOWING='PRECEDING' END IF C DO 134 J=JBGN,JEND,J3 C IF (KBND(J).NE.NB .OR. NSTR(J).LT.0) GO TO 134 C here if bands match. IF (DTYPE(J).EQ.'Y') THEN C this is sky; reset it. SKY=SIGNAL(J) PHNOISE=PHVAR(J) ERREST=ESTER(J) IF (KDIAPHRAGM.NE.-1) THEN C flag sky diaphragm code; expect star to match. DSKY=DIAPH(J) END IF ELSE C this is star; subtract sky. C IF (SKY.EQ.-3.E33) THEN C ERROR: no previous sky! CALL SPACE2 CARD=STARS(NSTR(J)) CARD(LWORD(CARD)+2:LWORD(CARD)+6)='at' WRITE(CARD(LWORD(CARD)+2:),'(F7.4)')DJOBS(J) CALL TV(CARD) CARD='!! Star not '//PRECEDED//' by sky !!' CALL TVN(CARD) CARD='Do you want to use the '//FOLLOWING//' sky?' CALL ASK(CARD,A) C IF (A.EQ.'Y' .OR. A.EQ.'O')THEN C look for following sky: DO 132 JJ=J+J3,JEND,J3 IF (KBND(JJ).NE.NB .OR. 1 NSTR(JJ).LT.0) GO TO 132 IF(DTYPE(JJ).EQ.'Y') THEN IF(KDIAPHRAGM.NE.-1 .AND. 1 DIAPH(JJ).NE.DIAPH(J)) GOTO 132 C fix this one, but do not set SKY. SIGNAL(J)=SIGNAL(J)-SIGNAL(JJ) C accumulate photon variances. PHVAR(J)=PHVAR(J)+PHVAR(JJ) IF (ESTER(J).NE.0.) 1 ESTER(J)=SQRT(ESTER(J)**2 + ESTER(JJ)**2) GO TO 134 END IF 132 CONTINUE CALL TV('No sky found! Quitting.') CALL STSEPI ELSE GO TO 999 END IF C ELSE C probably normal case. Check DIAPH: IF (KDIAPHRAGM.NE.-1 .AND. 1 DIAPH(J).NE.DSKY) THEN CALL SPACE2 CARD=STARS(NSTR(J)) CARD(LWORD(CARD)+2:LWORD(CARD)+6)='at' WRITE(CARD(LWORD(CARD)+2:), 1 '(F7.4)')DJOBS(J) CALL TV(CARD) CALL TV( 1 'Sky and Star taken through DIFFERENT diaphragms.') CALL TVN('Sky subtraction impossible.') GO TO 999 END IF C NORMAL case. Subtract sky: SIGNAL(J)=SIGNAL(J)-SKY PHVAR(J)=PHVAR(J)+PHNOISE IF(ESTER(J).NE.0.) ESTER(J)=SQRT(ESTER(J)**2 + ERREST**2) END IF C END IF C 134 CONTINUE C C Normal return. C RETURN C C C error return: C 999 CALL ASK('Do you want to try another sky-subtraction method?',A) IF (A.EQ.'Y' .OR. A.EQ.'O') THEN NERROR=1 RETURN ELSE IF (A.EQ.'H') THEN CALL TV('The photometric system changes with') CALL TVN('diaphragm size, so you cannot combine') CALL TVN('data taken through different holes.') CALL RTNCON('Try the "Nearest" option.',25) NERROR=1 RETURN ELSE CALL TV('Quitting.') CALL STSEPI END IF C END SUBROUTINE HALO(NB,DTYPE,KDIAPHRAGM,J,JBGN,JEND, SKYUSE) C C Picks the NEAREST STAR sample to sky(J), and returns it in SKYUSE. C C Arguments: C C NB is band number C DTYPE is 'Y' for sky, 'S' for star C KDIAPHRAGM is column pointer for diaphragm data C JBGN is observation number for start of night C JEND is observation number for end of night C C IMPLICIT NONE C C INTEGER NB, KDIAPHRAGM, JBGN, JEND C INCLUDE 'MID_REL_INCL:obs.inc' C CHARACTER DTYPE(MXOBS) C C C REAL SKY, SKY2, SKYAIR, SKYAZ, SKYTIM, PI,TWOPI, 1 SKYAIR2, SKYAZ2, SKYTIM2, S1, S2, SKYUSE C INTEGER J, JJ C C Meaning of local variables: C C Previous Following C -------- --------- C Star SKY SKY2 C Airmass SKYAIR SKYAIR2 C Azimuth SKYAZ SKYAZ2 C Time SKYTIM SKYTIM2 C Diaphragm DSKY DSKY2 C C INCLUDE 'MID_REL_INCL:mstars.inc' C PARAMETER (MSTARS=1650) C MSTARS is max.number of star-catalog entries. C commons for star catalog: CHARACTER *32 STARS COMMON /SCATA/ STARS(MSTARS) C CHARACTER DSKY*4, DSKY2*4 C DATA PI/3.1415926535/ C C C ***** BEGIN ***** C TWOPI=PI+PI SKY=-3.E33 SKY2=-3.E33 DSKY=' ' DSKY2=' ' SKYAIR2 = 0.0 SKYAZ2 = 0.0 SKYTIM2 = 0.0 C C J is current sky datum in band NB. C IF (DTYPE(J).EQ.'Y') THEN C this is sky. IF (KDIAPHRAGM.NE.-1) THEN C save new diaphragm. DSKY=DIAPH(J) ELSE C no diaphragm info. END IF ELSE CALL TV('Called HALO with non-sky datum.') RETURN END IF C C This is sky; look for PREVIOUS star. 20 IF (SKY.LT.0.) THEN C No previous star. Find one. DO 28 JJ=J-1,JBGN,-1 IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 28 IF(DTYPE(JJ).EQ.'S') THEN IF(KDIAPHRAGM.NE.-1 .AND. 1 DIAPH(JJ).NE.DIAPH(J)) GOTO 28 C OK; JJ is the previous star we seek. DSKY=DIAPH(JJ) SKY=SIGNAL(JJ) SKYAIR=AIRM(JJ) SKYAZ=AZIM(JJ) SKYTIM=DJOBS(JJ) GO TO 40 END IF 28 CONTINUE C No previous star. SKY=-3.E33 DSKY=' ' ELSE C probably normal case. Check DIAPH: IF (KDIAPHRAGM.NE.-1 .AND. DIAPH(J).NE.DSKY) THEN C Mismatch; force search for previous star. SKY=-3.E33 GO TO 20 ELSE C NORMAL case. Existing star OK. END IF END IF C C SKY is now valid previous value. Look for FOLLOWING value: C 40 IF (SKY2.LT.0.) THEN C Look for next star. DO 48 JJ=J+1,JEND IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 48 IF(DTYPE(JJ).EQ.'S') THEN IF(KDIAPHRAGM.NE.-1 .AND. 1 DIAPH(JJ).NE.DIAPH(J)) GOTO 48 C OK; JJ is the following star we seek. DSKY2=DIAPH(JJ) SKY2=SIGNAL(JJ) SKYAIR2=AIRM(JJ) SKYAZ2=AZIM(JJ) SKYTIM2=DJOBS(JJ) GO TO 60 END IF 48 CONTINUE C No following star. SKY2=-3.E33 DSKY2=' ' ELSE C existing sky2 might still be OK. Check it: IF (SKYTIM2.GT.DJOBS(J)) THEN C Time is OK. Check diaphragm: IF (KDIAPHRAGM.NE.-1) THEN IF (DIAPH(J).EQ.DSKY2) THEN C we are OK; go on. ELSE C need a new SKY2. SKY2=-3.E33 GO TO 40 END IF ELSE C No problem. Go on. END IF ELSE C We are past the old second star. Move it back & get new one: SKY=SKY2 DSKY=DSKY2 SKYAIR=SKYAIR2 SKYAZ=SKYAZ2 SKYTIM=SKYTIM2 SKY2=-3.E33 GO TO 40 END IF END IF C C Now compare to find "nearer" star: C 60 IF (SKY.GT.0. .AND. SKY2.GT.0.) THEN C Everything is hunky-dory. Find closer sample: S1=ABS(SKYTIM-DJOBS(J))*20. + ABS(SKYAIR-AIRM(J)) IF (ABS(SKYAZ-AZIM(J)).GT.PI) THEN S1=S1+ (SKYAIR+AIRM(J))*(TWOPI-ABS(SKYAZ-AZIM(J))) ELSE S1=S1+ (SKYAIR+AIRM(J))*ABS(SKYAZ-AZIM(J)) END IF S2=ABS(SKYTIM2-DJOBS(J))*20. + ABS(SKYAIR2-AIRM(J)) IF (ABS(SKYAZ2-AZIM(J)).GT.PI) THEN S2=S2+ (SKYAIR2+AIRM(J))*(TWOPI-ABS(SKYAZ2-AZIM(J))) ELSE S2=S2+ (SKYAIR2+AIRM(J))*ABS(SKYAZ2-AZIM(J)) END IF C IF (S2.LT.S1) THEN SKYUSE=SKY2 IF (S2.GT.0.1) THEN C complain. C CALL SKYFAR() END IF ELSE SKYUSE=SKY IF (S1.GT.0.1) THEN C complain. END IF END IF ELSE IF (SKY.GT.0.) THEN C SKY is the only good one. Use it. SKYUSE=SKY ELSE IF (SKY2.GT.0.) THEN C SKY2 is the only good one. Use it. SKYUSE=SKY2 ELSE C NEITHER star is any good. C set value to zero and return: SKYUSE=0. RETURN END IF C C OK, we have now rigged SKYUSE to be the value to use. Report it. C C C Normal return. C RETURN C END SUBROUTINE SKYNBR(NB,DTYPE,KDIAPHRAGM,JBGN,JEND) C C Picks the NEAREST sky sample, and subtracts it from star. C C Arguments: C C NB is band number C DTYPE is 'Y' for sky C KDIAPHRAGM is column pointer for diaphragm data C JBGN is observation number for start of night C JEND is observation number for end of night C C IMPLICIT NONE C C INTEGER NB, KDIAPHRAGM, JBGN, JEND C INCLUDE 'MID_REL_INCL:obs.inc' C CHARACTER DTYPE(MXOBS) C C C REAL SKY, SKY2, SKYAIR, SKYAZ, SKYTIM, PI,TWOPI, 1 PHNOISE, ERREST, SKYAIR2, SKYAZ2, SKYTIM2, 2 PHN2, ERR2, S1, S2, SKYUSE, PHNUSE, ERRUSE C INTEGER J, JJ C C Meaning of local variables: C C Previous Following C -------- --------- C Sky SKY SKY2 C Airmass SKYAIR SKYAIR2 C Azimuth SKYAZ SKYAZ2 C Time SKYTIM SKYTIM2 C Phot.noisePHNOISE PHN2 C Est.error ERREST ERR2 C Diaphragm DSKY DSKY2 C C INCLUDE 'MID_REL_INCL:mstars.inc' C PARAMETER (MSTARS=1650) C MSTARS is max.number of star-catalog entries. C commons for star catalog: CHARACTER *32 STARS COMMON /SCATA/ STARS(MSTARS) C CHARACTER CARD*48, DSKY*4, DSKY2*4 C DATA PI/3.1415926535/ C C TWOPI=PI+PI SKY=-3.E33 SKY2=-3.E33 DSKY=' ' DSKY2=' ' C DO 100 J=JBGN,JEND C IF (KBND(J).NE.NB .OR. NSTR(J).LT.0) GO TO 100 C here if bands match. C IF (DTYPE(J).EQ.'Y') THEN C this is sky. IF (KDIAPHRAGM.NE.-1) THEN C save new diaphragm. DSKY=DIAPH(J) ELSE C no diaphragm info. END IF SKY=SIGNAL(J) SKYAIR=AIRM(J) SKYAZ=AZIM(J) SKYTIM=DJOBS(J) PHNOISE=PHVAR(J) ERREST=ESTER(J) GO TO 100 ELSE C this is star; look for PREVIOUS sky. 20 IF (SKY.LT.0.) THEN C No previous sky. Find one. DO 28 JJ=J-1,JBGN,-1 IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 28 IF(DTYPE(JJ).EQ.'Y') THEN IF(KDIAPHRAGM.NE.-1 .AND. 1 DIAPH(JJ).NE.DIAPH(J)) GOTO 28 C OK; JJ is the previous sky we seek. DSKY=DIAPH(JJ) SKY=SIGNAL(JJ) SKYAIR=AIRM(JJ) SKYAZ=AZIM(JJ) SKYTIM=DJOBS(JJ) PHNOISE=PHVAR(J) ERREST=ESTER(J) GO TO 40 END IF 28 CONTINUE C No previous sky. SKY=-3.E33 DSKY=' ' ELSE C probably normal case. Check DIAPH: IF (KDIAPHRAGM.NE.-1 .AND. DIAPH(J).NE.DSKY) THEN C Mismatch; force search for previous sky. SKY=-3.E33 GO TO 20 ELSE C NORMAL case. Existing sky OK. END IF END IF END IF C C SKY is now valid previous value. Look for FOLLOWING value: C 40 IF (SKY2.LT.0.) THEN C Look for next sky. DO 48 JJ=J+1,JEND IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 48 IF(DTYPE(JJ).EQ.'Y') THEN IF(KDIAPHRAGM.NE.-1 .AND. 1 DIAPH(JJ).NE.DIAPH(J)) GOTO 48 C OK; JJ is the following sky we seek. DSKY2=DIAPH(JJ) SKY2=SIGNAL(JJ) SKYAIR2=AIRM(JJ) SKYAZ2=AZIM(JJ) SKYTIM2=DJOBS(JJ) PHN2=PHVAR(JJ) ERR2=ESTER(JJ) GO TO 60 END IF 48 CONTINUE C No following sky. SKY2=-3.E33 DSKY2=' ' ELSE C existing sky2 might still be OK. Check it: IF (SKYTIM2.GT.DJOBS(J)) THEN C Time is OK. Check diaphragm: IF (KDIAPHRAGM.NE.-1) THEN IF (DIAPH(J).EQ.DSKY2) THEN C we are OK; go on. ELSE C need a new SKY2. SKY2=-3.E33 GO TO 40 END IF ELSE C No problem. Go on. END IF ELSE C We are past the old second sky. Move it back & get new one: SKY=SKY2 DSKY=DSKY2 SKYAIR=SKYAIR2 SKYAZ=SKYAZ2 SKYTIM=SKYTIM2 PHNOISE=PHN2 ERREST=ERR2 SKY2=-3.E33 GO TO 40 END IF END IF C C Now compare to find "nearer" sky: C 60 IF (SKY.GT.0. .AND. SKY2.GT.0.) THEN C Everything is hunky-dory. Find closer sample: S1=ABS(SKYTIM-DJOBS(J))*20. + ABS(SKYAIR-AIRM(J)) IF (ABS(SKYAZ-AZIM(J)).GT.PI) THEN S1=S1+ (SKYAIR+AIRM(J))*(TWOPI-ABS(SKYAZ-AZIM(J))) ELSE S1=S1+ (SKYAIR+AIRM(J))*ABS(SKYAZ-AZIM(J)) END IF S2=ABS(SKYTIM2-DJOBS(J))*20. + ABS(SKYAIR2-AIRM(J)) IF (ABS(SKYAZ2-AZIM(J)).GT.PI) THEN S2=S2+ (SKYAIR2+AIRM(J))*(TWOPI-ABS(SKYAZ2-AZIM(J))) ELSE S2=S2+ (SKYAIR2+AIRM(J))*ABS(SKYAZ2-AZIM(J)) END IF C IF (S2.LT.S1) THEN SKYUSE=SKY2 PHNUSE=PHN2 ERRUSE=ERR2 IF (S2.GT.0.1) THEN C complain. C CALL SKYFAR() END IF ELSE SKYUSE=SKY PHNUSE=PHNOISE ERRUSE=ERREST IF (S1.GT.0.1) THEN C complain. END IF END IF ELSE IF (SKY.GT.0.) THEN C SKY is the only good one. Use it. SKYUSE=SKY PHNUSE=PHNOISE ERRUSE=ERREST ELSE IF (SKY2.GT.0.) THEN C SKY2 is the only good one. Use it. SKYUSE=SKY2 PHNUSE=PHN2 ERRUSE=ERR2 ELSE C AARRRGHHH!! NEITHER sky is any good! CARD='Can''t find a sky sample to match '//STARS(NSTR(J)) CALL TV(CARD) WRITE(CARD,'(''at '',F7.4)') DJOBS(J) IF (KDIAPHRAGM.NE.-1) CARD(30:)='in diaphragm '//DIAPH(J) CALL TVN(CARD) CALL TV('This value will be REJECTED.') CALL SPACE C set value to zero and continue: SKYUSE=SIGNAL(J) PHVAR(J)=1.E10 ERRUSE=1.E10 END IF C C OK, we have now rigged SKYUSE to be the value to subtract. Do it. C SIGNAL(J)=SIGNAL(J)-SKYUSE PHVAR(J)=PHVAR(J)+PHNUSE IF (ESTER(J).NE.0.) ESTER(J)=SQRT(ESTER(J)**2 +ERRUSE**2) C 100 CONTINUE C C Normal return. C RETURN C END SUBROUTINE REDSUB(BANDS,KMIN,KMAX,LASTK,MAXK,DTYPE) C C Subtracts red-leak data from corresponding bands. C Picks the NEAREST red-leak sample, and subtracts it from star. C C IMPLICIT NONE C REAL RLFACT, RL, RLVAR, ERVAR INTEGER KMIN, KMAX, LASTK, MAXK, NBRL, K, KRL, JBGN, NIGHT, 1 KK, J, JMIN, JMAX, JEND, ISTR, JOFF, J1, J2 C INCLUDE 'MID_REL_INCL:obs.inc' C INCLUDE 'MID_REL_INCL:mstars.inc' C PARAMETER (MSTARS=1650) C MSTARS is max.number of star-catalog entries. C commons for star catalog: CHARACTER *32 STARS COMMON /SCATA/ STARS(MSTARS) C CHARACTER DTYPE(MXOBS) C INCLUDE 'MID_REL_INCL:mbands.inc' C PARAMETER (MBANDS=9) INTEGER MXCOLR PARAMETER (MXCOLR=3*MBANDS) CHARACTER *8 BANDS(MXCOLR) C C Common for red leaks: C COMMON /REDLK/ NBRL(MXCOLR), RLFACT(MXCOLR) C C C Useful bands run from KMIN to LASTK. C Measured bands run from KMIN to KMAX, including red leaks. C MAXK marks end of expected red-leak bands, starting at NBANDS+1. C C C EXAMPLE of red-leak storage: C C In UBVRI system, NBANDS = 5. Suppose we measure only B,V,R, & BRL. C C K BAND(K) NBRL(K) C - ------- ------- C 1 U 0 C 2 B 6 <-- KMIN = 2 C 3 V 0 C 4 R 0 <-- LASTK = 4 C 5 I 0 C 6 BRL 0 <-- KMAX = MAXK = 6 C C C C DO 1000 K=KMIN,LASTK IF (NBRL(K).EQ.0) GO TO 1000 KRL=NBRL(K) JBGN=1 C DO 900 NIGHT=1,NIGHTS KK=0 DO 10 J=JBGN,NDATA C get limits for current night & passband. IF (NITE(J).NE.NIGHT) GO TO 20 IF (DTYPE(J).EQ.'Y') GO TO 10 KK=KK+1 IF (KK.EQ.1) JMIN=J JMAX=J 10 CONTINUE J=NDATA+1 20 JEND=J-1 C C JBGN, JEND index the first & last data of ANY kind for the night. C JMIN, JMAX index the first & last STAR data for the night & band. C KK counts star data for the night. C C DO 100 J=JMIN,JMAX C IF (DTYPE(J).EQ.'Y' .OR. NSTR(J).LT.0) GO TO 100 C here if datum is stellar. C IF (KBND(J).EQ.K) THEN C C this observation requires red-leak correction. ISTR=NSTR(J) C scan forward & back for nearest red-leak datum for same star: DO 50 JOFF=1,50 C look ahead... J1=J+JOFF IF (J1.GT.JMAX) GO TO 40 IF (KBND(J1).NE.KRL .OR. NSTR(J1).NE.ISTR) GO TO 40 RL=SIGNAL(J1) RLVAR=PHVAR(J1) ERVAR=ESTER(J1) GO TO 60 40 CONTINUE C look back... J2=J-JOFF IF (J2.LT.JMIN) GO TO 50 IF (KBND(J2).NE.KRL .OR. NSTR(J2).NE.ISTR) GO TO 50 RL=SIGNAL(J2) RLVAR=PHVAR(J2) ERVAR=ESTER(J2) GO TO 60 50 CONTINUE CALL TV('No red-leak datum found for') CALL TVN(STARS(ISTR)) C 60 CONTINUE IF (RL.LT.0.) THEN RL=0. CALL TV('Negative red-leak found for') CALL TVN(STARS(NSTR(J))) GO TO 100 END IF C C Correct for red leak: SIGNAL(J)=SIGNAL(J)-RL*RLFACT(K) PHVAR(J)=PHVAR(J)+RLVAR*RLFACT(K) IF(ESTER(J).NE.0.)ESTER(J)=SQRT(ESTER(J)**2+ERVAR*RLFACT(K)) C END IF C 100 CONTINUE C C prepare for next pass. JBGN=JEND+1 900 CONTINUE C 1000 CONTINUE C C Normal return. C RETURN C END SUBROUTINE NEED(LINES) C C Tell it how many lines you need for a plot. It reads the LOG C keyword and adds spaces to start a new page if necessary. C C IMPLICIT NONE C INTEGER LINES, NGOT, LINE, IUNIT, NULLS, ISTAT, NLPP, I C C CALL STKRDI('LOG',6,1, NGOT,LINE,IUNIT,NULLS,ISTAT) IF (ISTAT.NE.0) CALL TERROR(9999,'Could not read Keyword LOG') CALL STKRDI('LOG',7,1, NGOT,NLPP,IUNIT,NULLS,ISTAT) IF (ISTAT.NE.0) CALL TERROR(9999,'Could not read Keyword LOG') IF(LINE+LINES.GT.NLPP)THEN C we must space to a new page. DO 1 I=1,NLPP-LINE CALL TVN(' ') 1 CONTINUE ELSE C we are OK. END IF C RETURN END