C     PROGRAM RCVR1 EXAMPLE                                                     C                                                                               C     JOHN M. HOLT                                                              C     MIT LINCOLN LABORATORY                                                    C     MILLSTONE HILL FIELD STATION                                              C     P. O. BOX 73                                                              C     LEXINGTON, MASSACHUSETTS 02173                                            C                                                                               C     SAMPLE PROGRAM TO RECOVER INCOHERENT SCATTER DATA DESCRIBED IN            C     J. V. EVANS AND J. M. HOLT, MILLSTONE HILL THOMSON SCATTER                C     RESULTS FOR 1970, TECHNICAL REPORT 522, 11 MAY 1976. PREPARED             C     FOR THE NATIONAL SCIENCE FOUNDATION UNDER NSF GRANT GA-42230.             C                                                                               C     THIS PROGRAM ILLUSTRATES HOW SUBROUTINES CREAD AND RCVR1 MIGHT            C     BE USED TO RECOVER INCOHERENT SCATTER DATA FROM THE INSCON TAPE.          C     THE TAPE MUST INITIALLY BE AT THE BEGINNING OF THE THIRD FILE.            C     THE DAY IS SPECIFIED BY IDAY=700820 AND THE TYPE OF DATA IS               C     SPECIFIED BY KDT=101. THESE PARTICULAR VALUES SPECIFY THE                 C     ELECTRON DENSITY DATA COLLECTED ON 23-24 MARCH, 1970. BOTH OF             C     THESE NUMBERS ARE FROM THE INSCON SUMMARY, WHICH LISTS THE                C     POLYNOMIAL FITS APPEARING ON THE INSCON TAPE. THE TIME AND                C     ALTITUDE ARE SPECIFIED BY TIME=3600.0*(12.0+24.0)-RELTME AND              C     Z=405.0. THESE CORRESPOND TO 1200 EST ON 24 MARCH AND 405.0 KM.           C     NOTE THAT THE SPECIFIED TIME AND ALTITUDE NEED NOT BE THE TIME            C     AND ALTITUDE OF A MEASUREMENT. THEY NEED ONLY FALL WITHIN THE             C     TIME AND ALTITUDE RANGE OF THE DATA. IN THE EXPRESSION FOR TIME,          C     3600.0 IS A FACTOR CONVERTING FROM HOURS TO SECONDS, 12.0 IS THE          C     HOUR (EST) AND 24.0 INDICATES THE SECOND DAY OF THE RUN. RELTME           C     IS ON THE INSCON FIT TAPE AND CONVERTS TO THE TIME ORIGIN USED            C     IN THE POLYNOMIAL FITS.                                                   C                                                                               C     THE DAY ON WHICH THE EXPERIMENT BEGINS, 23 MARCH (DAY 82) IN THIS         C     CASE, IS ASSUMED TO BEGIN AT 0000 GMT. SOME CARE MUST THEREFORE           C     BE TAKEN WHEN AN EXPERIMENT BEGINS LESS THAN 5 HOURS BEFORE               C     MIDNIGHT EST. IF THE EXPERIMENT BEGAN AT 2130 EST ON                      C     16 JULY (EST), THE FIRST DAY LISTED ON THE INSCON SUMMARY                 C     WOULD BE 17 JULY, SINCE THE EXPERIMENT BEGAN AT 0230 GMT ON               C     17 JULY (GMT). IN THAT CASE, RELTME=75600.=3600.*21.                      C                                                                               C     THE 23-24,JUNE,1970 EXPERIMENT BEGAN JUST BEFORE 2400 GMT ON              C     23 JUNE. THE TIME GIVEN IN THE INSCON SUMMARY, 1901 EST, IS               C     SLIGHTLY AFTER 0000 GMT ON 24 JUNE BECAUSE THE LISTED TIME                C     CORRESPONDS TO THE BEGINNING OF THE FIRST C-MODE OF THE                   C     EXPERIMENT SEQUENCE. THE FIRST A-MODE OF THIS EXPERIMENT BEGAN            C     AT 2347 GMT ON 23 JUNE.                                                   C                                                                                     COMMON/CDAT/ IDAYNO, KDAT, NX, NY, NPX, NPY, XT, YT, DMAX, DMIN,               1             XMAX, XMIN, YMAX, YMIN, RELTME, SHIFT, GAMX, DELX,                2             FACX, GAMY, DELY, FACY, A, SIG, ISTAT                              DIMENSION XT(64), YT(64), GAMX(32), DELX(32), FACX(32), GAMY(32),              1          DELY(32), FACY(32), A(32,32), SIG(64), ISTAT(64,64)             C                                                                               C     .....SPECIFY DAY, DATA TYPE, TIME AND ALTITUDE.....                             IDAY = 700820                                                                   KDT = 101                                                                 C                                                                               C     .....READ DATA SETS UNTIL DESIRED FIT IS ENCOUNTERED.....                    10 CALL CREAD                                                                      IF(IDAYNO .EQ. IDAY .AND. KDAT .EQ. KDT) GO TO 20                               GO TO 10                                                                  C                                                                               C     .....RECOVER DATA AND DERIVATIVES AT SPECIFIED TIME AND ALTITUDE..           20 T = 12.0                                                                        TIME = 3600.*(T+24.) - RELTME                                                   Z = 405.                                                                        CALL RCVR1(TIME, Z, P, PDX, PDY)                                          C                                                                               C     .....OUTPUT RESULTS.....                                                        WRITE(6,101) T, Z, P, PDX, PDY                                              101 FORMAT(7H TIME =, F6.2/ 11H ALTITUDE =, F7.1/ 4H P =, E12.5/                   1       8H DP/DT =, E12.5/ 8H DP/DZ =, E12.5)                                    STOP                                                                      C                                                                                     END                                                                                                                                                                                                                                       C     PROGRAM TEST RCVR1                                                TEST0010C                                                                       TEST0020C     INSCON IS A FORTRAN PROGRAM WHICH CALCULATES LEAST MEAN SQUARE    TEST0030C     FITS TO, AND PLOTS CONTOUR DIAGRAMS OF, INCOHERENT SCATTER DATA   TEST0040C     OBTAINED AT MILLSTONE HILL. THE FIT VALUE AND FIRST DERIVATIVE    TEST0050C     VALUES CORRESPONDING TO A GIVEN TIME AND ALTITUDE MAY BE RECOVEREDTEST0060C     BY MEANS OF SUBROUTINES CREAD AND RCVR1. FOR EXAMPLE THIS PROGRAM TEST0070C     FIRST CALLS CREAD TO READ THE CARD DECK GENERATED BY INSCON. IT   TEST0080C     THEN CALLS RCVR1 TO RECOVER THE INSCON FIT VALUES AND DERIVATIVES TEST0090C     AT THE EXPERIMENTAL TIMES AND ALTITUDES.                          TEST0100C                                                                       TEST0110      COMMON/CDAT/ IDAYNO, KDAT, NX, NY, NPX, NPY, XT, YT, DMAX, DMIN,  TEST0120     1             XMAX, XMIN, YMAX, YMIN, RELTME, SHIFT, GAMX, DELX,   TEST0130     2             FACX, GAMY, DELY, FACY, A, SIG, ISTAT                TEST0140      DIMENSION XT(64), YT(64), GAMX(32), DELX(32), FACX(32), GAMY(32), TEST0150     1          DELY(32), FACY(32), A(32,32), SIG(64), ISTAT(64,64)     TEST0160C                                                                       TEST0170      CALL CREAD                                                        TEST0180      WRITE(6,101)                                                      TEST0190      DO 10 I=1,NPX                                                     TEST0200      ITIME = XT(I) + RELTME                                            TEST0210      IF(ITIME .GT. 86400) ITIME=ITIME-86400                            TEST0220      IHRS = ITIME/3600                                                 TEST0230      MINS = (ITIME-IHRS*3600)/60                                       TEST0240      ITIME = 100*IHRS + MINS                                           TEST0250      DO 10 J=1,NPY                                                     TEST0260      CALL RCVR1(XT(I), YT(J), P, PDX, PDY)                             TEST0270      WRITE(6,201) ITIME, YT(J), P, PDX, PDY                            TEST0280   10 CONTINUE                                                          TEST0290      STOP                                                              TEST0300C                                                                       TEST0310C     .....FORMAT STATEMENTS.....                                       TEST0320  101 FORMAT(1H1, 4X, 4HTIME, 5X, 8HALTITUDE, 8X, 1HP, 13X, 5HDP/DX,    TEST0330     1       11X, 5HDP/DY)                                              TEST0340  201 FORMAT(5X, I4, 5X, F6.0, F16.5, 4X, E12.5, 4X, E12.5)             TEST0350C                                                                       TEST0360      END                                                               TEST0370                                                                                      SUBROUTINE RCVR1(X, Y, P, PDX, PDY)                               RCV10010C                                                                       RCV10020C     JOHN M. HOLT                                                      RCV10030C     MIT LINCOLN LABORATORY                                            RCV10040C     MILLSTONE HILL FIELD STATION                                      RCV10050C     P. O. BOX 73                                                      RCV10060C     LEXINGTON, MASSACHUSETTS 02173                                    RCV10070C                                                                       RCV10080C     RCVR1 USES THE INSCON FIT PARAMETERS STORED IN CDAT TO RECOVER    RCV10090C     THE INSCON FIT VALUE P , AND ITS DERIVATIVES WITH RESPECT TO X ANDRCV10100C     Y, PDX AND PDY, AT TIME X AND ALTITUDE Y. KDAT INDICATES          RCV10110C     THE TYPE OF DATA, FOR EXAMPLE ELECTRON DENSITY. THIS INFORMATION  RCV10120C     IS ALSO ON THE FIRST CARD OF THE INSCON OUTPUT DECK. X IS IN      RCV10130C     SECONDS, WITH X=0 AT THE HOUR PRECEDING THE FIRST EXPERIMENTAL    RCV10140C     DATA(RELTME). THE SMALLEST AND LARGEST PERMISSIBLE VALUES OF      RCV10150C     X ARE GIVEN BY XMIN AND XMAX. Y IS IN KILOMETERS AND MAY RANGE    RCV10160C     FROM YMIN TO YMAX. IF X OR Y IS OUTSIDE THE ALLOWED RANGE, RCVR1  RCV10170C     RETURNS P=0. THE NPX EXPERIMENTAL VALUES OF X ARE STORED IN XT.   RCV10180C     THE NPY EXPERIMENTAL VALUES OF Y ARE STORED IN YT. THE DATA TO    RCV10190C     WHICH INSCON FITS AN NX BY NY LEAST MEAN SQUARE POLYNOMIAL RANGE  RCV10200C     FROM DMIN TO DMAX. EACH TERM OF THE POLYNOMIAL IS OF THE FORM     RCV10210C     PX(I)*PY(J)*A(J,I). THE A(J,I) ARE COEFFICIENTS OUTPUT BY INSCON. RCV10220C     THE PX(I) AND PY(J) ARE POLYNOMIALS WHICH ARE ORTHONORMAL OVER    RCV10230C     XT(I) AND YT(J) RESPECTIVELY. THEY ARE DEFINED BY RECURSION       RCV10240C     COEFFICIENTS GAMX, DELX, FACX, GAMY, DELY, FACY. IF DESIRED A     RCV10250C     LOWER ORDER LEAST MEAN SQUARE FIT MAY BE RECOVERED BY DECREASING  RCV10260C     NX AND/OR NY. ISTAT(J,I) IS AN NPY BY NPX ARRAY WHICH INDICATES   RCV10270C     THE QUALITY OF THE FIT AT (XT(I), YT(J)). IF ISTAT(J,I)=1,        RCV10280C     THE FIT IS GOOD. IF ISTAT(J,I)=0, THE FIT MAY BE BAD AND RCVR1    RCV10290C     RETURNS P=-P. LARGE TIME GAPS IN THE DATA ARE INDICATED BY        RCV10300C     ISTAT(J,I)=2. RCVR1 RETURNS P=-P IF X IS WITHIN THE GAP.          RCV10310C     SIG(J) IS THE SQUARE ROOT OF THE MEAN SQUARE DEVIATION ABOUT      RCV10320C     THE FIT OF THE DATA AT ALTITUDE YT(J). SHIFT IS A POSITIVE        RCV10330C     CONSTANT WHICH IS ADDED BEFORE FITTING TO DATA WHICH TAKES ON     RCV10340C     NEGATIVE VALUES.                                                  RCV10350C                                                                       RCV10360      COMMON/CDAT/ IDAYNO, KDAT, NX, NY, NPX, NPY, XT, YT, DMAX, DMIN,  RCV10370     1             XMAX, XMIN, YMAX, YMIN, RELTME, SHIFT, GAMX, DELX,   RCV10380     2             FACX, GAMY, DELY, FACY, A, SIG, ISTAT                RCV10390      DIMENSION XT(64), YT(64), GAMX(32), DELX(32), FACX(32), GAMY(32), RCV10400     1          DELY(32), FACY(32), A(32,32), SIG(64), ISTAT(64,64)     RCV10410      DIMENSION PX(32), PY(32), PX1(32), PY1(32)                        RCV10420C                                                                       RCV10430      IF(X .LE. XMAX .AND. X .GE. XMIN .AND.                            RCV10440     1   Y .LE. YMAX .AND. Y .GE. YMIN) GO TO 10                        RCV10450      P = 0.                                                            RCV10460      RETURN                                                            RCV10470   10 XP = (2.*X-XMAX-XMIN)/(XMAX-XMIN)                                 RCV10480      YP = (2.*Y-YMAX-YMIN)/(YMAX-YMIN)                                 RCV10490      CALL SET1(XP, NX, GAMX, DELX, FACX, PX, PX1)                      RCV10500      CALL SET1(YP, NY, GAMY, DELY, FACY, PY, PY1)                      RCV10510      P = 0.                                                            RCV10520      PDX = 0.                                                          RCV10530      PDY = 0.                                                          RCV10540      DO 20 I=1,NX                                                      RCV10550      DO 20 J=1,NY                                                      RCV10560      P = P + PX(I)*PY(J)*A(J,I)                                        RCV10570      PDX = PDX + PX1(I)*PY(J)*A(J,I)                                   RCV10580   20 PDY = PDY + PX(I)*PY1(J)*A(J,I)                                   RCV10590      P = ((DMAX-DMIN)*P+DMAX+DMIN)/2.                                  RCV10600      PDX = (DMAX-DMIN)*PDX/(XMAX-XMIN)                                 RCV10610      PDY = (DMAX-DMIN)*PDY/(YMAX-YMIN)                                 RCV10620      IF(IPL(X,Y) .EQ. 3) P=-P                                          RCV10630      RETURN                                                            RCV10640C                                                                       RCV10650      END                                                               RCV10660                                                                                      SUBROUTINE CREAD                                                  CRED0010C                                                                       CRED0020C     CREAD READS THE CARD DECK GENERATED BY INSCON. THE CONTENTS OF    CRED0030C     THE CARDS ARE STORED IN LABELED COMMON BLOCK CDAT. INSCON PUNCHES CRED0040C     ISTAT(J,I) ONLY IF, GIVEN J, EITHER THERE IS AN I SUCH THAT       CRED0050C     ISTAT(J,I) IS NOT 0, OR J=NPY. IF NO CARD IS PRESENT FOR A GIVEN  CRED0060C     J, CREAD SETS ISTAT(J,I)=1.                                       CRED0070C                                                                       CRED0080      COMMON/CDAT/ IDAYNO, KDAT, NX, NY, NPX, NPY, XT, YT, DMAX, DMIN,  CRED0090     1             XMAX, XMIN, YMAX, YMIN, RELTME, SHIFT, GAMX, DELX,   CRED0100     2             FACX, GAMY, DELY, FACY, A, SIG, ISTAT                CRED0110      DIMENSION XT(64), YT(64), GAMX(32), DELX(32), FACX(32), GAMY(32), CRED0120     1          DELY(32), FACY(32), A(32,32), SIG(64), ISTAT(64,64),    CRED0130     2          IST(64)                                                 CRED0140C                                                                       CRED0150      READ(5,101) IDAYNO, KDAT, NX, NY, NPX, NPY                        CRED0160      READ(5,201) (XT(I), I=1,NPX)                                      CRED0170      READ(5,201) (YT(J), J=1,NPY)                                      CRED0180      READ(5,201) DMAX, DMIN, XMAX, XMIN, YMAX, YMIN, RELTME, SHIFT     CRED0190      N = MAX0(NX,NY)                                                   CRED0200      READ(5,201) (GAMX(I), DELX(I), FACX(I),                           CRED0210     1             GAMY(I), DELY(I), FACY(I), I=1,N)                    CRED0220      DO 10 I=1,NX                                                      CRED0230   10 READ(5,201) (A(J,I), J=1,NY)                                      CRED0240      READ(5,201) (SIG(J), J=1,NPY)                                     CRED0250      J1 = 1                                                            CRED0260   20 READ(5,301) YP, (IST(I), I=1,NPX)                                 CRED0270      DO 60 J=J1,NPY                                                    CRED0280      IF(ABS(YT(J)-YP) .LT. .1 .OR. YT(J) .GE. 9999.9) GO TO 40         CRED0290      DO 30 I=1,NPX                                                     CRED0300   30 ISTAT(J,I) = 1                                                    CRED0310      GO TO 60                                                          CRED0320   40 DO 50 I=1,NPX                                                     CRED0330   50 ISTAT(J,I) = IST(I)                                               CRED0340      IF(J .EQ. NPY) RETURN                                             CRED0350      J1 = J + 1                                                        CRED0360      GO TO 20                                                          CRED0370   60 CONTINUE                                                          CRED0380C                                                                       CRED0390  101 FORMAT(73X, I6/ 5(1X, I3))                                        CRED0400  201 FORMAT(6(1X, E12.5))                                              CRED0410  301 FORMAT(F7.1, 4X, 64I1)                                            CRED0420C                                                                       CRED0430      END                                                               CRED0440                                                                                      SUBROUTINE SET1(X, N, GAM, DEL, FAC, P, P1)                       SET10010C                                                                       SET10020C     SET1 CALCULATES THE FIRST N ORTHONORMAL POLYNOMIALS               SET10030C     P (DEFINED BY GAM, DEL, FAC) AND THEIR DERIVATIVES P1 AT          SET10040C     ARGUMENT X.                                                       SET10050C                                                                       SET10060      DIMENSION P(1), P1(1), GAM(1), DEL(1), FAC(1)                     SET10070      P(1) = 1.                                                         SET10080      P(2) = X - GAM(2)                                                 SET10090      P1(1) = 0.                                                        SET10100      P1(2) = 1.                                                        SET10110      DO 10 L=3,N                                                       SET10120      P(L) = (X-GAM(L))*P(L-1) - DEL(L)*P(L-2)                          SET10130   10 P1(L) = P(L-1) + (X-GAM(L))*P1(L-1) - DEL(L)*P1(L-2)              SET10140      DO 20 L=1,N                                                       SET10150      P(L) = P(L)/FAC(L)                                                SET10160   20 P1(L) = P1(L)/FAC(L)                                              SET10170      RETURN                                                            SET10180      END                                                               SET10190                                                                                      FUNCTION IPL(X,Y)                                                 IPL 0010C                                                                       IPL 0020C     IPL DETERMINES WHETHER THE INSCON FIT IS GOOD AT POINT (X,Y).     IPL 0030C                                                                       IPL 0040      COMMON/CDAT/ IDAYNO, KDAT, NX, NY, NPX, NPY, XT, YT, DMAX, DMIN,  IPL 0050     1             XMAX, XMIN, YMAX, YMIN, RELTME, SHIFT, GAMX, DELX,   IPL 0060     2             FACX, GAMY, DELY, FACY, A, SIG, ISTAT                IPL 0070      DIMENSION XT(64), YT(64), GAMX(32), DELX(32), FACX(32), GAMY(32), IPL 0080     1          DELY(32), FACY(32), A(32,32), SIG(64), ISTAT(64,64)     IPL 0090C                                                                       IPL 0100      DX = ABS(X-XT(1))                                                 IPL 0110      DO 10 I=2,NPX                                                     IPL 0120      DX1 = ABS(X-XT(I))                                                IPL 0130      IF(DX1 .GT. DX) GO TO 20                                          IPL 0140   10 DX = DX1                                                          IPL 0150   20 IP = I - 1                                                        IPL 0160      DY = ABS(Y-YT(1))                                                 IPL 0170      DO 30 J=2,NPY                                                     IPL 0180      DY1 = ABS(Y-YT(J))                                                IPL 0190      IF(DY1 .GT. DY) GO TO 40                                          IPL 0200   30 DY = DY1                                                          IPL 0210   40 JP = J - 1                                                        IPL 0220      IF(ISTAT(JP,IP)-1) 50,60,70                                       IPL 0230   50 IPL = 3                                                           IPL 0240      RETURN                                                            IPL 0250   60 IPL = 2                                                           IPL 0260      RETURN                                                            IPL 0270   70 IPL = 2                                                           IPL 0280      IF(X .GT. XT(IP) .AND. ISTAT(JP,IP+1) .EQ. 2 .OR.                 IPL 0290     1   X .LT. XT(IP) .AND. ISTAT(JP,IP-1) .EQ. 2) IPL=3               IPL 0300      RETURN                                                            IPL 0310      END                                                               IPL 0320