C *********** C **PREPARE** C *********** UNB00030 C VERSION JULY 1975 C PREPARE MJV DATA FOR INPUT TO ADJUSTMENT UNB00040 C PROCESS PASS BY PASS AS FOLLOWS - UNB00050 C READ MJV DATA (MJVIN) UNB00060 C SCALE KEPLER PARAMETERS (SKEP) UNB00070 C SCALE EPHEMERAL PARAMETERS (SEPH) UNB00080 C APPLY IONOSPHERIC CORRECTION TO DOPPLERS (IONOS) UNB00090 C COMPUTE PASS ALERT (ALERT) UNB00100 C FIT BASE FUNCTIONS TO EPHEMERAL PARAMETERS (FITEPH) UNB00110 C COMPUTE SATELLITE XYZ AT 4.6 SECOND INTERVALS (SATXYZ) UNB00120 C APPLY TROPOSPHERIC CORRECTION TO DOPPLERS (TROPO) UNB00130 C CONVERT SAT COORDS TO TCA RANGE/VELOCITY COORD SYSTEM (RANVEL) UNB00140 C COMPUTE 2D FIX ON GUIER PLANE (GUIER) UNB00150 C CARD INPUTS DESCRIBED IN SUBR CARD1,CARD3,CARD5,CARD12,CARD15 C MAXIMUM OF 18 CARDS REQUIRED C DATA VECTORS C DOP1 = HIGH CHANNEL DOPPLERS C DOP2 = LOW CHANNEL DOPPLERS C DOP3 = TIME C DOP4 = VACUUM COUNTS C DOP5 = IONOSPHERIC CORRECTION C DOP6 = TROPO CORRECTION C DOP7 = ELEVATION ANGLE IN DEGREES C DOP8 = EDIT STRING C DOP9 = RESIDUALS C**MACHINE DEPENDENT STATEMENT UNB00170 IMPLICIT REAL*8(A-H,O-Z) UNB00180 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00190 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,7),DOP8(260),DOP9(260) DIMENSION PASSK(16),KPASS(16),NULLA(20),XS(3) DIMENSION IDOP(20) EQUIVALENCE (FHEAD(54),PCODE),(PHEAD(1),PASS) $ ,(FHEAD(86),TAU1N),(FHEAD(113),TAU2N),(FHEAD(185),PASSK(1)) $ ,(PHEAD(26),DOPNUM) C READ INPUT CARDS CALL HEADER DO 10 I = 1,16 10 KPASS(I) = PASSK(I) DO 15 I = 1,20 IDOP(I) = 0 15 NULLA(I) = 0 C INITIALIZE PRINTOUT HEADING IF(ILIST .EQ. 0) WRITE(IPR,2008) IF(KPASS(1) .EQ. 0) GO TO 60 UNB00250 C PROGRAM MODE = PROCESS UP TO 16 SELECTED PASSES J2 = 0 DO 30 I = 1,16 IF(KPASS(I) .EQ. 0) GO TO 40 ILIST = 3 J1 = J2 + 1 J2 = KPASS(I) DO 20 J = J1,J2 NPASS = NPASS + 1 PASS = NPASS NULL = 0 CALL PINIT CALL MJVIN IF(IEDS .NE. 0) GO TO 80 20 CONTINUE CALL SKEP CALL SEPH IF(NULL .NE. 0) GO TO 30 CALL FITEPH IF(NULL .NE. 0) GO TO 30 CALL TVECT(1) CALL DOPTCA IF(NULL .NE. 0) GO TO 30 CALL ALERT IF(NULL .NE. 0) GO TO 30 IF(PCODE .EQ. 2.) CALL POLE CALL METOBS CALL RANVEL IF(TAU1N .NE. TAU2N) CALL COMBIN ITC = FHEAD(60) CALL TVECT(4 + ITC) CALL EDIT IF(NULL .NE. 0) GO TO 30 CALL IONOS CALL TROPO CALL GUIER IF(NULL .NE. 0) GO TO 30 CALL SPASS 30 CONTINUE 40 CONTINUE RETURN C PROGRAM MODE = PROCESS ALL PASSES IN DATA SET 60 NPASS = NPASS + 1 PASS = NPASS NULL = 0 CALL PINIT CALL MJVIN IF(IEDS .NE. 0) GO TO 80 CALL SKEP CALL SEPH IF(NULL .NE. 0) GO TO 70 CALL FITEPH IF(NULL .NE. 0) GO TO 70 CALL TVECT(1) CALL DOPTCA IF(NULL .NE. 0) GO TO 70 CALL ALERT IF(NULL .NE. 0) GO TO 70 IF(PCODE .EQ. 2.) CALL POLE CALL METOBS CALL RANVEL IF(TAU1N .NE. TAU2N) CALL COMBIN ITC = FHEAD(60) CALL TVECT(4 + ITC) CALL EDIT IF(NULL .NE. 0) GO TO 70 CALL IONOS CALL TROPO CALL GUIER IF(NULL .NE. 0) GO TO 70 CALL SPASS NDOP = DOPNUM DO 65 I = 1,NDOP J = DOP8(I) + 11. 65 IDOP(J) = IDOP(J) + 1 70 NULLA(NULL+1) = NULLA(NULL+1) + 1 GO TO 60 80 CONTINUE WRITE(IPR,1001) NPASS = NPASS - 1 WRITE(IPR,1002) NPASS,(NULLA(I),I=1,7) WRITE(IPR,1003) IDOP(12),(IDOP(I),I=1,10) RETURN UNB00810 C FORMAT STATEMENTS 1001 FORMAT(1X,15HEND OF DATA ) 1002 FORMAT(1H1,4X,I5,13H TOTAL PASSES//5X,I5,9H ACCEPTED//5X,I5, $ 11H MISSED TCA/5X,I5,13H TOO FEW DOPS/5X,I5,8H BAD FIX/5X,I5, $ 10H BAD ALERT/5X,I5,10H BAD ORBIT/5X,I5,11H BAD LOCKON) 1003 FORMAT(//5X,I10,18H DOPPLERS ACCEPTED//5X,I10,11H MISCLOSURE/5X, $ I10,9H RESIDUAL/5X,I10,7H COMMON/5X,I10,9H SYMMETRY/5X,I10, $ 17H TOO FAR FROM TCA/5X,I10,6H SLOPE/5X,I10,10H ELEVATION/5X,I10 $ ,11H LO CHANNEL/5X,I10,11H HI CHANNEL/5X,I10,13H LOSS OF LOCK/) 2008 FORMAT(10H1PASS KEPT,20X,9HSAT EL AZ,9X,21HINJECTION LOCK ON , UNB02140 * 8HDOPPLERS/) UNB02150 END UNB00820 SUBROUTINE HEADER C READ INPUT DATA CARDS AND COMPUTE ELEMENTS OF ARRAY FHEAD C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) EQUIVALENCE (FHEAD(54),PCODE),(PHEAD(1),PASS) CALL CARDA CALL CARDB CALL CARDC CALL CARDD IF(PCODE .EQ. 2.) CALL CARDE CALL SHEAD RETURN END SUBROUTINE CARDA C INITIALIZE DATA SET REFERENCE NUMBERS, CONSTANTS AND POINTERS C AND READ INPUT CARDS ONE AND TWO C**MACHINE DEPENDENT STATEMENT UNB00040 IMPLICIT REAL*8(A-H,O-Z) UNB00050 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00060 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) COMMON /FORM/ IFORM DIMENSION TITLE(20) C DATA SET REFERENCE NUMBERS UNB00250 C MET=WEATHER INPUT, IN=MJV INPUT, IOUT=OUTPUT RESULTS, IOUT2=SCRATCH,UNB00260 C ICR=CARD READER, IPR=PRINTER1, IPN=PUNCH, IPR2=PRINTER2, IPR3=PRINT3UNB00270 MET = 1 UNB00280 IN = 2 UNB00290 IN = 5 IOUT = 3 UNB00300 IOUT2= 4 UNB00310 ICR = 5 UNB00320 IPR = 6 UNB00330 IPN = 7 UNB00340 IPR2 = 8 UNB00350 IPR3 = 9 UNB00360 C MACHINE DEPENDENT PARAMETERS. VEL=VACUUM VELOCITY OF PROP (M/SEC) UNB00370 C WE=ROTATION RATE OF EARTH (RAD/MIN). FO=GROUND REFERENCE FREQ (HZ) PI = 3.141592653589793D0 UNB00390 VEL = 299792.5D3 UNB00400 WE = 4.3752695D-3 UNB00410 FO = 400D6 C DERIVED CONSTANTS UNB00430 RADG = PI / 180D0 UNB00440 WAVENO = FO / VEL C BIT PERIODS FOR BR(SATELLITE) AND CBR(RECEIVER) TIMING BITS = 120D0 / 6103D0 BITR = 7865000D0 / FO C POINTERS UNB00460 IKEPT = 0 UNB00490 IEDS = 0 UNB00500 C RESET FILE HEADER ARRAY DO 10 I = 1,200 10 FHEAD(I) = 0. CARD ONE CONTAINS LISTING OPTION AND INPUT FORMAT(2I1) C ILIST = LEVEL OF OUTPUT LISTING UNB00780 C 0 - ONE LINE SUMMARY OF EACH PASS ONLY UNB00790 C 1 - SUMMARY PLUS ERROR MESSAGES UNB00800 C 2 - COMPLETE LISTING OF COMPUTED DATA UNB00810 C 3 - COMPLETE INPUT AND OUTPUT LISTING UNB00820 C IFORM = 1 FOR PROGRAM MAJORITY OUTPUT C 2 FOR GEODETIC SURVEY OF CANADA MAJORITY VOTE OUTPUT C ADD CODE TO SUBR MJVIN FOR ADDITIONAL FORMATS READ(ICR,1002) ILIST,IFORM C SET UP DEFAULT VALUES IF(IFORM .EQ. 0) IFORM = 1 IF(ILIST .GE. 1) WRITE(IPR,1004) CARD TWO DESCRIBES THIS COMPUTER RUN IN UP TO 80 CHARACTERS (20A4) READ(ICR,1001) TITLE UNB00640 WRITE(IPR,2001) TITLE IF(IFORM .EQ. 1) WRITE(IPR,2002) IF(IFORM .EQ. 2) WRITE(IPR,2003) RETURN C FORMAT STATEMENTS 1001 FORMAT(20A4) UNB01890 1002 FORMAT(2I1) 1004 FORMAT(1X,30HSUBROUTINE CARDA ENTERED ) 2001 FORMAT(1H1,50X,30HPROGRAM **PREPARE** OUTPUT /20X,20A4) 2002 FORMAT(5X,50HFORMAT OF INPUT DATA = PROGRAM MAJORITY OUTPUT ) 2003 FORMAT(5X,50HFORMAT OF INPUT DATA = GEODETIC SURVEY MJV FORMAT ) END SUBROUTINE CARDB C READ INPUT CARDS THREE AND FOUR C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION COORD(3),CSIG(3),X0(3),R0(3),XG(3),PLH(3),SXYZ(3,3), $ SPLH(3,3),TOPROT(3,3) DIMENSION R(3),NAX(3) DATA NAX(1),NAX(2),NAX(3)/3,2,-2/ EQUIVALENCE (FHEAD(1),STN),(FHEAD(2),XG(1)),(FHEAD(5),SXYZ(1,1)), $ (FHEAD(14),PLH(1)),(FHEAD(17),SPLH(1,1)),(FHEAD(26),TOPROT(1,1)) $ ,(FHEAD(35),RG),(FHEAD(36),A),(FHEAD(37),B),(FHEAD(38),X0(1)), $ (FHEAD(41),R0(1)),(FHEAD(44),S0) C**MACHINE DEPENDENT STATEMENT UNB00180 SQRT(PI) = DSQRT(PI) UNB00210 IF(ILIST .GE. 1) WRITE(IPR,1000) CARD THREE DESCRIBES STATION (A4,4X,8F8.0) C STNAM = UP TO 4 CHAR STATION NAME C STN = STATION NUMBER C CCODE = COORDINATE CODE C 1 IF COORD(3) AND CSIG(3) REFER TO XYZ COORDS (IN M) C 2 IF COORD(3) AND CSIG(3) REFER TO LAT/LONG/HT COORDS C (LAT/LONG IN DEGREES + DECIMAL DEGREES, HT IN M) C COORD(3) = STATION COORDINATES C CSIG(3) = STANDARD DEVIATIONS OF STATION COORDINATES (IN M) READ(ICR,1002) STNAM,STN,CCODE,COORD,CSIG C SET UP DEFAULTS IF(STN .EQ. 0.) STN = 1. IF(CCODE .EQ. 0.) CCODE = 1. DO 5 I = 1,3 5 IF(CSIG(I) .EQ. 0.) CSIG(I) = 100. C CHECK FOR NONZERO COORDINATES IF(COORD(1).EQ.0. .OR. COORD(2).EQ.0. .OR. (CCODE .NE.2. .AND. $ COORD(3).EQ.0.)) GO TO 40 WRITE(IPR,2002) STNAM,STN CARD FOUR DESCRIBES DATUM (9F8.0) C A = SEMI-MAJOR AXIS (IN M) C B = SEMI-MINOR AXIS (IN M) OR RECIPROCAL OF FLATTENING C X0(3) = TRANSLATION COMPONENTS (IN M) C RO(3) = ROTATION COMPONENTS (IN ARCSEC) C S0 = SCALE DIFFERENCE (IN PPM) READ(ICR,1003) A,B,X0,R0,S0 C SET UP DEFAULTS IF(A .EQ. 0.) A = 6378388D0 IF(B .EQ. 0.) B = 297D0 WRITE(IPR,2003) A,B,X0,R0,S0 C CONVERT COORDINATES AND PRINT IN BOTH FORMS ICOORD = CCODE GO TO (10,20),ICOORD 10 DO 15 I = 1,3 UNB01490 SXYZ(I,I) = CSIG(I)**2 15 XG(I) = COORD(I) UNB01500 CALL RECGEO(A,B,X0,R0,S0,XG,SXYZ,PLH,SPLH) GO TO 30 UNB01480 C IF INPUT IS LAT/LONG/HT CONVERT AND STORE AS XYZ 20 DO 25 I = 1,3 SPLH(I,I) = CSIG(I)**2 25 PLH(I) = COORD(I) IF(PLH(2) .LT. 0.) PLH(2) = PLH(2) + 360. CALL GEOREC(A,B,X0,R0,S0,PLH,SPLH,XG,SXYZ) 30 CONTINUE RG = SQRT(XG(1)**2 + XG(2)**2 + XG(3)**2) UNB01720 WRITE(IPR,2006) XG,PLH,RG,((SXYZ(I,J),J=1,3),(SPLH(I,J),J=1,3), $ I=1,3) C COMPUTE TOPROT = TRANSFORMATION MATRIX FROM AVG TERR TO TOPOCENTRIC PHI = PLH(1) * RADG XLAMB = PLH(2) * RADG R(1) = XLAMB - PI R(2) = PHI - PI / 2. R(3) = 0. CALL ROTREF(3,NAX,R,TOPROT) RETURN C ERROR RETURN 40 WRITE(IPR,2011) STOP C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE CARDB ENTERED ) 1002 FORMAT(A4,4X,8F8.0) 1003 FORMAT(10F8.0) 2002 FORMAT(20X,17HOBSERVING STATION ,5X,A4,5X,10HNUMBER ,F10.0/) 2003 FORMAT(5X,25HDATUM PARAMETERS A ,F20.3/29X,1HB,F20.3/28X, $ 2HX0,F20.3/28X,2HY0,F20.3/28X,2HZ0,F20.3/,28X,2HR1,F20.3/28X, $ 2HR2,F20.3/28X,2HR3,F20.3/28X,2HDS,F20.3/) 2006 FORMAT(5X,14HAPPROX COORDS ,3F13.3,5X,2F13.5,F13.3,5X,F13.3// $ 3(19X,3F13.3,5X,3F13.3/)) 2011 FORMAT(1X,35HZERO STATION COORDINATES DETECTED ) END SUBROUTINE CARDC C READ CARDS FIVE TO ELEVEN C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION SEPHSD(3),VEPHSD(3),PASS(16) EQUIVALENCE (FHEAD(45),ECODE),(FHEAD(46),SEPHSD(1)),(FHEAD(49), $ VEPHSD(1)),(FHEAD(52),POLY),(FHEAD(53),APPSIG),(FHEAD(54),PCODE) $ ,(FHEAD(55),YEAR),(FHEAD(56),RCODE),(FHEAD(57),SERIAL),(FHEAD(58 $ ),DELAY),(FHEAD(59),DELSIG),(FHEAD(60),TCODE),(FHEAD(61),DOPSIG) $ ,(FHEAD(62),FGO),(FHEAD(63),DFG),(FHEAD(64),DAYFGO),(FHEAD(65), $ FGOSIG),(FHEAD(66),DFGSIG),(FHEAD(67),ELCUT),(FHEAD(68),CUT), $ (FHEAD(69),SLOPEM),(FHEAD(70),YMISCL),(FHEAD(71),ALF1),(FHEAD(72 $ ),ELMIN),(FHEAD(73),DOPMIN),(FHEAD(74),XITER),(FHEAD(75),EPS1), $ (FHEAD(76),EPS2),(FHEAD(77),EPS3),(FHEAD(78),EPS4),(FHEAD(79), $ EPS5),(FHEAD(80),ALF2),(FHEAD(81),WCODE),(FHEAD(82),TK), $ (FHEAD(83),PSL),(FHEAD(84),PWV),(FHEAD(85),TROPSD) $ ,(FHEAD(185),PASS(1)) IF(ILIST .GE. 1) WRITE(IPR,1000) CARD FIVE DESCRIBES EPHEMERIS(10F8.0) C ECODE = EPHEMERIS CODE C 1 FOR PRECISE EPHEMERIS C 2 FOR BROADCAST EPHEMERIS C SEPHSD(3) = EPHEMERIS SHIFT STANDARD DEVIATIONS C VEPHSD(3) = EPHEMERIS VELOCITY STANDARD DEVIATIONS C POLY = NUMBER OF FUNCTIONS USED TO APPROXIMATE EPHEMERALS C APPSIG = MAXIMUM STD DEV FOR EPHEMERAL APPROXIMATION (IN BR. UNITS) C PCODE = POLAR MOTION CODE C 1 EPHEMERIS REFERRED TO CIO. NO POLE COORDS READ, NO CORRECTION C 2 EPHEMERIS REFERRED TO INST POLE. POLE COORDINATES READ IN READ(ICR,1003) ECODE,SEPHSD,VEPHSD,POLY,APPSIG,PCODE C SET UP DEFAULTS IF(ECODE .EQ. 0.) ECODE = 2. IF(SEPHSD(1) .EQ. 0.) SEPHSD(1) = 2. IF(SEPHSD(2) .EQ. 0.) SEPHSD(2) = 0.5 IF(SEPHSD(3) .EQ. 0.) SEPHSD(3) = 1. IF(POLY .EQ. 0.) POLY = 4. IF(APPSIG .EQ. 0.) APPSIG = 0.25 IF(PCODE .EQ. 0.) PCODE = 1. IF(ECODE .EQ. 1.) WRITE(IPR,2007) SEPHSD,VEPHSD IF(ECODE .EQ. 2.) WRITE(IPR,2008) SEPHSD,VEPHSD WRITE(IPR,2018) POLY,APPSIG IF(PCODE .EQ. 1.) WRITE(IPR,2011) IF(PCODE .EQ. 2.) WRITE(IPR,2012) CARD SIX DESCRIBES MEASUREMENT INSTRUMENT (7F8.0) C YEAR = YEAR IN WHICH MEASUREMENTS WERE MADE C RCODE = RECEIVER CODE C 1 = MARCONI 722 WITH 7 DIGIT :OPPLERS C 2 = ITT 5001 (UNB CONFIGURATION) C 3 = MAGNAVOX 702CA C 4 = MARCONI 722B WITH 8 DIGIT DOPPLERS (ONE FRACTIONAL DIGIT) C 5 = MARCONI 722B WITH 9 DIGIT DOPPLERS (TWO FRACTIONAL DIGITS) C SERIAL = RECEIVER SERIAL NUMBER C DELAY = RECEIVER NOMINAL DELAY (IN MICROSECONDS) C DELSIG = DELAY STANDARD DEVIATION (IN MICROSECONDS) C TCODE = MEASUREMENT TIMING CODE C 1 = SATELLITE TIME BASE (BR DATA) C 2 = RECEIVER TIME BASE (CBR DATA) C DOPSIG = DOPPLER MEASUREMENT STANDARD DEVIATION (IN COUNTS) READ(ICR,1003) YEAR,RCODE,SERIAL,DELAY,DELSIG,TCODE,DOPSIG C SET UP DEFAULTS IF(YEAR .EQ. 0.) YEAR = 1975. IF(RCODE .EQ. 0.) RCODE = 5. IF(DELAY .EQ. 0.) DELAY = 550. IF(DELSIG .EQ. 0.) DELSIG = 100. IF(TCODE .EQ. 0.) TCODE = 2. IF(DOPSIG .EQ. 0.) DOPSIG = 2. WRITE(IPR,2004) YEAR,RCODE,SERIAL,DELAY,DELSIG,DOPSIG IF(TCODE .EQ. 1.) WRITE(IPR,2009) IF(TCODE .EQ. 2.) WRITE(IPR,2010) CARD SEVEN DESCRIBES RECEIVER OSCILLATOR (5F8.0) C MODEL IS FG(T) = FGO + DFG * (T - DAYFGO) WHERE C T = TIME OF PASS IN DAYS C FGO = FREQ OFFSET ON DAYFGO (IN PARTS IN 10(10)) C DFG = FREQ DRIFT (IN PARTS IN 10(10) PER DAY) C DAYFGO = DATUM DAY FOR FGO C FGOSIG = FREQ OFFSET STANDARD DEVIATION (IN PARTS IN 10(10)) C DFGSIG = FREQ DRIFT STANDARD DEVIATION (IN PARTS IN 10(10) /DAY) READ(ICR,1003) FGO,DFG,DAYFGO,FGOSIG,DFGSIG C SET UP DEFAULTS IF(FGOSIG .EQ. 0.) FGOSIG = 1000. IF(DFGSIG .EQ. 0.) DFGSIG = 10. WRITE(IPR,2005) FGO,DFG,DAYFGO,FGOSIG,DFGSIG CARD EIGHT CONTAINS DOPPLER REJECTION OPTIONS (5F8.0) C ELCUT = ELEVATION CUTOFF ANGLE IN DEG (DEFAULT = 7.5) C CUT = MAXIMUM NUMBER OF INTERVALS FROM TCA (DEFAULT = 200.) C SLOPEM = MINIMUM SLOPE OF DOPPLER CURVE IN HZ/S (DEFAULT = 0.) C YMISCL = MAXIMUM MISCLOSURE VALUE IN COUNTS (DEFAULT = 1000.) C ALF1 = SIGNIF LEVEL FOR NORMALITY TEST (DEFAULT = 99 PERCENT) READ(ICR,1003) ELCUT,CUT,SLOPEM,YMISCL,ALF1 C SET UP DEFAULT VALUES IF(ELCUT .EQ. 0.) ELCUT = 7.5 IF(CUT .EQ. 0.) CUT = 200. IF(YMISCL .EQ. 0.) YMISCL = 1000. IF(ALF1 .EQ. 0.) ALF1 = 0.99 WRITE(IPR,2014) ELCUT,CUT,SLOPEM,YMISCL,ALF1 CARD NINE CONTAINS PASS REJECTION OPTIONS (9F8.0) C ELMIN = MINIMUM PASS ELEVATION IN DEGREES (DEFAULT = 15.) C DOPMIN = MINIMUM NUMBER OF DOPPLERS PER PASS(DEFAULT =4.) C XITER = MAXIMUM NUMBER OF ITERATIONS TO CONVERGE (DEFAULT = 5) C EPS1 = CONVERGENCE TEST FOR COORDINATES (DEFAULT = 1E-3) C EPS2 = CONVERGENCE TEST FOR FREQUENCY OFFSET (DEFAULT = 0.1) C EPS3 = MAXIMUM CORRECTION TO A PRIORI COORDINATES (DEFAULT=1000.) C EPS4 = MAXIMUM STD DEV OF COMPUTED COORDINATES (DEFAULT = 100.) C EPS5 = MAXIMUM STD DEV OF COMPUTED FREQ OFFSET (DEFAULT = 10.) C ALF2 = SIGNIF LEVEL FOR VARIANCE FACTOR CHI-SQ TEST (DEFAULT=99PC) READ(ICR,1003) ELMIN,DOPMIN,XITER,EPS1,EPS2,EPS3,EPS4,EPS5,ALF2 C SET UP DEFAULT VALUES IF(ELMIN .EQ. 0.) ELMIN = 15. IF(DOPMIN .EQ. 0.) DOPMIN = 4. IF(XITER .EQ. 0.) XITER = 5. IF(EPS1 .EQ. 0.) EPS1 = 1E-3 IF(EPS2 .EQ. 0.) EPS2 = 0.1 IF(EPS3 .EQ. 0.) EPS3 = 1000. IF(EPS4 .EQ. 0.) EPS4 = 100. IF(EPS5 .EQ. 0.) EPS5 = 10. IF(ALF2 .EQ. 0.) ALF2 = 0.99 WRITE(IPR,2015) ELMIN,DOPMIN,XITER,EPS1,EPS2,EPS3,EPS4,EPS5,ALF2 CARD TEN DESCRIBES TROPO REFRACTION CORRECTION (5F8.0) C WCODE = SURFACE WEATHER CODE C 1 IF STANDARD VALUES ARE TO BE USED C 2 IF OBSERVED MET VALUES TO BE READ FROM TAPE C TK = STANDARD SURFACE TEMPERATURE (IN K) C PSL = STANDARD SEA LEVEL PRESSURE (IN MB) C PWV = STANDARD WATER VAPOUR PRESSURE (IN MB) C TROPOSD = TROPO SCALING STD DEV READ(ICR,1003) WCODE,TK,PSL,PWV,TROPSD C SET UP DEFAULT VALUES IF(WCODE .EQ. 0.) WCODE = 1. IF(TK .EQ. 0.) TK = 293. IF(PSL .EQ. 0.) PSL = 1014. IF(PWV .EQ. 0.) PWV = 20. IF(TROPSD .EQ. 0.) TROPSD = 0.1 IF(WCODE .EQ. 1.) WRITE(IPR,2016) IF(WCODE .EQ. 2.) WRITE(IPR,2017) WRITE(IPR,2006) TK,PSL,PWV,TROPSD CARD ELEVEN CONTAINS PASS SELECTION ARRAY (16F5.0) C KPASS = ARRAY SELECTING UP TO 16 PASSES TO BE PROCESSED C IF KPASS(1) = 0 ALL PASSES ARE PROCESSED UNB00960 READ(ICR,1004) PASS IF(PASS(1) .NE. 0.) WRITE(IPR,2013) PASS RETURN UNB01840 C FORMAT STATEMENTS UNB01880 1000 FORMAT(1X,30HSUBROUTINE CARDC ENTERED ) 1003 FORMAT(10F8.0) 1004 FORMAT(16F5.0) 2004 FORMAT(5X,6HYEAR ,F10.0/5X,14HRECEIVER CODE ,F10.0,10H SERIAL NU, $ 5HMBER ,F10.0/5X,14HNOMINAL DELAY ,F10.0,9H STD DEV ,F10.0/5X, $ 16HDOPPLER STD DEV ,F15.5,7H COUNTS) 2005 FORMAT(5X,41HRECEIVER FREQ OFFSET IN PARTS IN 10(10) = ,F10.3, UNB02090 $ 3H + ,F10.3,19H * (TIME IN DAYS - ,F11.3,2H )/37X,9HSTD DEVS , $ F10.3,3X,F10.3/) 2006 FORMAT(5X,14HSTD TEMP ,F10.2/5X,14HSTD SL-PRESS ,F10.2/ $ 5X,14HSTD WV-PRESS ,F10.2/5X,20HTROPO SCALING STDDEV ,F5.2) 2007 FORMAT(5X,24HPRECISE EPHEMERIS USED,2(/5X,8HSTD DEVS,3F15.5)) 2008 FORMAT(5X,24HBROADCAST EPHEMERIS USED,2(/5X,8HSTD DEVS,3F15.5)) 2009 FORMAT(5X,35HSATELLITE TIME BASE (BR) USED /) 2010 FORMAT(5X,35HRECEIVER TIME BASE (CBR) USED /) 2011 FORMAT(5X,50HNO POLAR MOTION CORRECTION MADE TO SAT COORDS /) 2012 FORMAT(5X,50HPOLE CORRECTION APPLIED TO SATELLITE COORDINATES /) 2013 FORMAT(5X,50HONLY THE FOLLOWING PASSES WILL BE PROCESSED $ /5X,16F5.0) 2014 FORMAT( 5X,13HREJECT DOPPLE, $ 5HRS LT,F5.1,4H DEG/21X,2HGT,F5.0,13H INT FROM TCA/21X, 4HIF S, $ 7HLOPE LT,F5.0,5H HZ/S/21X,11HIF MISCL GT,F10.1,7H COUNTS/21X, $ 8HIF FAILS,F5.2,17H % NORMALITY TEST/) 2015 FORMAT( 5X,18HREJECT PASS IF EL , $ 2HLT,F5.1/17X,5HIF LT,F5.0,5H DOPS/17X,5HIF GT,F5.0,5H ITER, $ 16H FOR CONVERGE TO ,2E10.3/17X,13HIF NOT WITHIN ,F10.3, BUG***** $ 18H OF APRIORI COORDS/17X,14HIF STD DEVS GT ,2F10.3/17X, $ 8HIF FAILS,F5.2,18H % CHI-SQUARE TEST/) 2016 FORMAT(5X,25HSTANDARD MET VALUES USED ) 2017 FORMAT(5X,25HOBSERVED MET VALUES USED ) 2018 FORMAT( 3X,F5.0,21H-TERM APPROX TO EPHEM /5X,16HREJECT APPROX IF, $ 33H CHI-SQUARE FAILS WITH APRIORI VF,F10.5) END UNB02200 SUBROUTINE CARDD C READ INPUT CARDS TWELVE TO FOURTEEN C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION TAU1(26),TAU2(26),TAUN(15) EQUIVALENCE (FHEAD(86),TAU1N),(FHEAD(87),TAU1(1)),(FHEAD(113), $ TAU2N),(FHEAD(114),TAU2(1)),(FHEAD(140),TAUN(1)) IF(ILIST .GE. 1) WRITE(IPR,1000) CARDS TWELVE AND THIRTEEN CONTAIN INTEGRATION INTERVALS OF INPUT DATA C INTERVALS EXPRESSED IN BIT PERIODS (APPROX 20 MS). C SET OF INTERVALS COMPRISING TWO MINUTES DATA ONLY, IS SPECIFIED. READ (ICR,1001) TAU1 C SET UP DEFAULTS IF(TAU1(1) .NE. 0.) GO TO 5 TAU1(1) = 1638. TAU1(2) = 1404. TAU1(3) = 1638. TAU1(4) = 1423. 5 CONTINUE C FIND NUMBER OF INTERVALS/TWO MINS, AND CHECK TOTAL BIT PERIODS=6103 T6103 = 0. DO 10 I = 1,26 IF(TAU1(I) .EQ. 0.) GO TO 15 10 T6103 = T6103 + TAU1(I) I = 27 15 TAU1N = I - 1 IF(T6103 .NE. 6103.) GO TO 100 IF(TAU1N .GT. 26.) GO TO 100 NTAU1 = TAU1N WRITE(IPR,1002) TAU1N,(TAU1(I),I=1,NTAU1) CARD FOURTEEN CONTAINS NUMBER OF INPUT INTEGRATION INTERVALS TO BE C COMBINED FOR EACH PROCESSING INTERVAL (16F5.0) C TAUN(1) .EQ. 0. TO PROCESS INPUT INTERVALS C TAUN(1) .GE. TAU1N TO PROCESS TWO MINUTE INTERVALS READ(ICR,1001) TAUN IF(TAUN(1) .EQ. 0.) GO TO 50 IF(TAUN(1) .GE. TAU1N) GO TO 60 C FIND NUMBER OF PROCESSING INTERVALS / TWO MINUTES, AND CHECK C TOTAL OF TAUN ARRAY EQUALS TAU1N TAUM = 0. DO 30 I = 1,15 IF(TAUN(I) .EQ. 0.) GO TO 35 30 TAUM = TAUM + TAUN(I) I = 16 35 TAU2N = I - 1 IF(TAUM .NE. TAU1N) GO TO 110 C COMPUTE TAU2 VECTOR K2 = 0 NTAU2 = TAU2N DO 40 I = 1,NTAU2 TAU2(I) = 0. K1 = K2 + 1 K2 = K1 + TAUN(I) - 1 DO 40 K = K1,K2 40 TAU2(I) = TAU2(I) + TAU1(K) GO TO 70 C PROCESS USING INPUT INTEGRATION INTERVALS 50 TAU2N = TAU1N NTAU2 = TAU2N DO 55 I = 1,NTAU2 55 TAU2(I) = TAU1(I) GO TO 70 C PROCESS USING TWO MINUTE INTEGRATION INTERVALS 60 TAU2N = 1. NTAU2 = TAU2N TAU2(1) = 6103. TAUN(1) = TAU1N 70 WRITE(IPR,1004) TAU2N,(TAU2(I),I=1,NTAU2) RETURN C ERROR RETURNS 100 WRITE(IPR,1005) T6103,TAU1N,TAU1 STOP 110 WRITE(IPR,1006) TAUM,TAUN STOP C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE CARDD ENTERED ) 1001 FORMAT(16F5.0) 1002 FORMAT(5X,F4.0,30H INPUT INTEGRATION INTERVALS ,2(/5X,13F6.0)) 1004 FORMAT(5X,F4.0,30H PROCESS INTEGRATION INTERVALS,2(/5X,13F6.0)) 1005 FORMAT(5X,35HINPUT INTEGRATION INTERVAL ERROR ,2F10.0,2(/5X, $ 13F6.0)) 1006 FORMAT(5X,40HPROCESSING INTEGRATION INTERVAL ERROR ,F10.0/5X, $ 15F4.0) END SUBROUTINE CARDE C READ INPUT CARDS FIFTEEN TO SEVENTEEN C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION DPOLE(10),XP(10),YP(10) EQUIVALENCE (FHEAD(155),DPOLE(1)),(FHEAD(165),XP(1)),(FHEAD(175), $ YP(1)) IF(ILIST .GE. 1) WRITE(IPR,1000) CARD FIFTEEN CONTAINS DAYNUMBERS FOR POLE COORDINATES CARD SIXTEEN CONTAINS X-COORDINATES OF POLE RELATIVE TO CIO CARD SEVENTEEN CONTAINS Y-COORDINATES OF POLE RELATIVE TO CIO READ(ICR,1007) DPOLE,XP,YP WRITE(IPR,2007) (DPOLE(I),XP(I),YP(I),I=1,10) RETURN C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE CARDE ENTERED ) 1007 FORMAT(10F5.0/10F8.2/10F8.2) 2007 FORMAT(5X,20H DAY XPOLE YPOLE/10(5X,F5.0,2F8.2/)) END SUBROUTINE GEOREC(A,B,X0,R0,S0,PLH,SPLH,XG,SXYZ) C GEODETIC TO AVERAGE TERRESTRIAL CONVERSION C INPUT GEODETIC COORDINATES ARE LAT/LONG/HT IN DEGREES AND M C REFERRED TO ELLIPSOID WITH AXES A,B C OUTPUT AVERAGE TERRESTRIAL COORDINATES ARE CARTESIAN, IN M C SEVEN PARAMETERS TRANSFORMATION PERFORMED C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS DIMENSION X(3),XG(3),X0(3),R0(3),PLH(3),SPLH(3,3),SXYZ(3,3), $ R(3),ROT(3,3),NAX1(3),NAX2(3) DATA NAX1(1),NAX1(2),NAX1(3)/1,2,3/ DATA NAX2(1),NAX2(2),NAX2(3)/-2,2,3/ SQRT(PI) = DSQRT(PI) SIN(PI) = DSIN(PI) COS(PI) = DCOS(PI) IF(ILIST .GE. 1) WRITE(IPR,1000) PHI = PLH(1) * RADG XLAMB = PLH(2) * RADG HT = PLH(3) CP = COS(PHI) SP = SIN(PHI) CL = COS(XLAMB) SL = SIN(XLAMB) BOA = (B/A)**2 IF(B .LT. 6D6) BOA = (1. - 1. / B) **2 XN = A / SQRT(CP**2 + BOA * SP**2) C GEODETIC TO CARTESIAN CONVERSION XG(1) = (XN + HT) * CP * CL XG(2) = (XN + HT) * CP * SL XG(3) = (XN * BOA + HT) * SP C THREE ROTATIONS DO 10 I = 1,3 R(I) = R0(I) * RADG / 3600. X(I) = XG(I) 10 XG(I) = 0D0 CALL ROTREF(3,NAX1,R,ROT) DO 20 I = 1,3 DO 20 J = 1,3 20 XG(I) = XG(I) + ROT(I,J) * X(J) C SCALE CHANGE DO 30 I = 1,3 30 XG(I) = XG(I) * (1D0 + S0*1D-6) C THREE TRANSLATIONS DO 40 I = 1,3 40 XG(I) = XG(I) + X0(I) C TRANSFORM COVARIANCE MATRIX R(1) = 0. R(2) = PI / 2. - PHI R(3) = PI - XLAMB CALL ROTREF(3,NAX2,R,ROT) CALL PROP(SXYZ,ROT,SPLH,3,3) RETURN C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE GEOREC ENTERED ) END SUBROUTINE RECGEO(A,B,X0,R0,S0,XG,SXYZ,PLH,SPLH) C AVERAGE TERRESTRIAL TO GEODETIC CONVERSION C INPUT AVERAGE TERRESTRIAL CARTESIAN COORDINATES IN M C SEVEN PARAMETER TRANSFORMATION PERFORMED C OUTPUT GEODETIC COORDINATES LAT/LONG/HT IN DEGREES AND M C REFERRED TO ELLIPSOID WITH AXES A,B C ITERATIVE COMPUTATION CONTINUED UNTIL CONVERGENCE TO 1E-10 RADIANS IN C PHI AND A*1E-10 METRES IN HT, OR FOR MAXIMUM OF 20 ITERATIONS C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS DIMENSION X(3),XX(3),XG(3),X0(3),R0(3),PLH(3),SPLH(3,3),SXYZ(3,3), $ R(3),ROT(3,3),NAX1(3),NAX2(3) DATA NAX1(1),NAX1(2),NAX1(3)/3,2,1/ DATA NAX2(1),NAX2(2),NAX2(3)/3,2,-2/ ABS(PI) = DABS(PI) SQRT(PI) = DSQRT(PI) SIN(PI) = DSIN(PI) COS(PI) = DCOS(PI) ATAN(PI) = DATAN(PI) ATAN2(X1,Y1) = DATAN2(X1,Y1) IF(ILIST .GE. 1) WRITE(IPR,1000) C INITIALIZE EPS2 = 1E-10 EPS1 = A * EPS2 ITER = 20 BOA = (B/A)**2 IF(B .LT. 6D6) BOA = (1. - 1. / B) **2 E2=1.-BOA C THREE TRANSLATIONS DO 10 I = 1,3 10 X(I) = XG(I) - X0(I) C SCALE CHANGE DO 20 I = 1,3 20 X(I) = X(I) / (1D0 + S0 * 1D-6) C THREE ROTATIONS DO 30 I = 1,3 30 R(I) = R0(I) * RADG / 3600. CALL ROTREF(3,NAX1,R,ROT) DO 40 I = 1,3 XX(I) = 0D0 DO 40 J = 1,3 40 XX(I) = XX(I) + ROT(I,J) * X(J) X1 = XX(1) Y1 = XX(2) Z1 = XX(3) C COMPUTE LONGITUDE NON-ITERATIVELY P = SQRT(X1**2 + Y1**2) XLAMB = ATAN2(Y1,X1) C INITIAL APPROXIMATIONS FOR ITERATION XN=A HT = SQRT(X1**2 + Y1**2 + Z1**2) - A PHI = ATAN((Z1/P)/(1.-E2*XN/(XN+HT))) C ITERATION LOOP DO 50 I = 1,ITER HT1=HT PHI1=PHI XN = A / SQRT(COS(PHI)**2 + BOA * SIN(PHI)**2) HT = P / COS(PHI) - XN PHI = ATAN((Z1/P)/(1.-E2*XN/(XN+HT))) IF(ABS(HT1-HT) .LT. EPS1 .AND. ABS(PHI1-PHI) .LT. EPS2) GO TO 60 50 CONTINUE GO TO 70 60 PLH(1) = PHI / RADG PLH(2) = XLAMB / RADG PLH(3) = HT C TRANSFORM COVARIANCE MATRIX R(1) = XLAMB - PI R(2) = PHI - PI / 2. R(3) = 0. CALL ROTREF(3,NAX2,R,ROT) CALL PROP(SPLH,ROT,SXYZ,3,3) RETURN C ERROR RETURN 70 WRITE(IPR,1001) STOP C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE RECGEO ENTERED ) 1001 FORMAT(1X,30HTOO MANY ITERATIONS IN RECGEO ) END SUBROUTINE PROP(QX,T,QY,NX,NY) C PROPAGATION OF COVARIANCE MATRICES C IF X = T Y C THEN QX = T QY T(T) C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) DIMENSION QX(NX,NX),T(NX,NY),QY(NY,NY) DO 10 I = 1,NX DO 10 J = 1,NX QX(I,J) = 0D0 DO 10 K = 1,NY DO 10 L = 1,NY 10 QX(I,J) = QX(I,J) + T(I,K) * T(J,L) * QY(K,L) RETURN END SUBROUTINE SHEAD C STORE AND PRINT CONTENTS OF ARRAY FHEAD IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) IF(ILIST .GE. 1) WRITE(IPR,1000) WRITE(IPN,1002) FHEAD IF(ILIST .GE. 3) WRITE(IPR,1002) FHEAD RETURN C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE SHEAD ENTERED ) 1002 FORMAT(1X,F6.0,3F13.3,3F10.0/1X,6F10.0,F13.8/1X,F13.8,F13.3,5F10.0 $ /1X,4F10.0,3F13.10/1X,6F13.10/1X,F13.3,2F10.1/1X,3F9.3,4F9.4/ $ 1X,F3.0,6F5.2,F3.0,F5.2,F3.0,F6.0,F3.0,F5.0,F6.0,F6.0,F3.0,F6.2/ $ 1X,F8.1,F6.1,F5.0,F8.0,F6.0,F5.1,F5.0, F5.1,F7.1,F5.2/ $ 1X,F5.1,F4.0,F4.0,F5.3,F4.1,F6.0,F5.0,F4.0,F5.2,F3.0,F5.0,F6.0, $ F4.0,F5.2/1X,F4.0,12F6.0/1X,13F6.0/1X,F6.0,F4.0,11F6.0/1X,13F6.0 $ /1X,2F6.0,15F4.0/1X,10F5.0,5F5.1/1X,15F5.1/1X,15F5.0/1X,F5.0) END SUBROUTINE PINIT C INITIALIZE PASS ARRAYS IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) IF(ILIST .GT. 0) WRITE(IPR,1000) DO 10 I = 1,100 10 PHEAD(I) = 0. DO 20 I = 1,20 ORB1(I) = 0. 20 ORB3(I) = 0. DO 30 I = 1,70 30 ORB2(I) = 0. DO 50 I = 1,260 DO 50 J = 1,9 50 DUMMY(I,J) = 0D0 RETURN 1000 FORMAT(1X,30HSUBROUTINE PINIT ENTERED ) END SUBROUTINE MJVIN C READ MAJORITY VOTED DATA IN PROGRAM **MAJORITY** OUTPUT FORMAT C IRCVR = INTEGER IDENTIFYING RECEIVER MAKE C NROW2 = NUMBER OF ROWS PER TWO-MINUTE MESSAGE C NDPMSG = NUMBER OF MESSAGES CONTAINING DOPPLER DATA (UP TO 10) C NVAR = NUMBER OF VARIABLE PARAMETER SETS MAJORITY VOTED (UP TO 17) C LOCK = INDEX OF VARIABLE PARAMETER SET REFERRING TO LOCK ON TIME C IFIX(14) = MAJORITY VOTED FIXED PARAMETERS C IVAR(NVAR) = MAJORITY VOTED VARIABLE PARAMETERS C IDOP(2,NROW2,NDPMSG) = HIGH AND LOW CHANNEL 4.6 SECOND DOPPLERS C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), $ DOP7(260),DOP8(260),DOP9(260) COMMON /FORM/ IFORM DIMENSION STAT(6),ELEM(6),IFIX(14),IVAR(17),IDOP(26,10,2) DIMENSION ERR(10) DATA R,S/1HR,1HS/ EQUIVALENCE (FHEAD(56),RCODE),(PHEAD(17),STAT(1)),(PHEAD(23), $ ELEM(1)),(PHEAD(51),XLOCK) $ ,(FHEAD(86),TAU1N) $ ,(FHEAD(57),SERIAL),(FHEAD(60),TCODE),(PHEAD(2),SAT),(PHEAD(3), $ DLOCK),(PHEAD(4),TLOCK) IF(ILIST .GE. 1) WRITE(IPR,1000) GO TO (5,100,150),IFORM 5 CONTINUE NTAU1 = TAU1N JRCVR = RCODE READ(IN,1001,END=50,ERR=50) IRCVR,NROW2,NDPMSG,NVAR,LOCK,IFIX, $ (IVAR(I),I=1,NVAR) IF(ILIST .GE. 1) WRITE(IPR,1003) IRCVR,NROW2,NDPMSG,NVAR,LOCK,IFIX $ ,(IVAR(I),I=1,NVAR) IF(IRCVR .NE. JRCVR) GO TO 60 IF(NROW2 .NE. NTAU1) GO TO 70 XLOCK = LOCK DO 10 I = 1,14 10 ORB1(I) = IFIX(I) DO 20 I = 1,NVAR 20 ORB2(I) = IVAR(I) STAT(1) = 1. STAT(2) = 2. ELEM(1) = 14. ELEM(2) = NVAR READ(IN,1002,END=50,ERR=50) (((IDOP(J,K,L),L=1,2),J=1,NROW2), $ K=1,NDPMSG) IF(ILIST .GE. 1) WRITE(IPR,1004) (((IDOP(J,K,L),L=1,2),J=1,NROW2), $ K = 1,NDPMSG) DO 30 L = 1,2 DO 30 K = 1,NDPMSG DO 30 J = 1,NROW2 N = J + (K-1)*NROW2 IF(L .EQ. 1) DOP1(N) = IDOP(J,K,L) / 100D0 IF(L .EQ. 2) DOP2(N) = IDOP(J,K,L) / 100D0 30 CONTINUE C COMPUTE NDOP NDOP = NROW2 * NDPMSG DO 35 I = 1,NDOP J = 1 + NDOP - I IF(DOP1(J) .NE. 0. .OR. DOP2(J) .NE. 0.) GO TO 37 35 CONTINUE 37 NDOP = J IF(ILIST .GE. 3) WRITE(IPR,1007) (I,DOP1(I),DOP2(I),I=1,NDOP) STAT(4) = 1. STAT(5) = 2. ELEM(4) = NDOP ELEM(5) = NDOP RETURN C END OF DATA SET 50 IEDS = 1 RETURN C ERROR RETURNS 60 WRITE(IPR,1005) JRCVR,IRCVR STOP 70 WRITE(IPR,1006) NTAU1,NROW2 STOP C READ DATA IN GEODETIC SURVEY OF CANADA MAJORITY VOTED FORMAT 100 CONTINUE NVAR = 14 XLOCK = 3. C READ HEADER, VARIABLE PARAMETERS, AND FIXED PARAMETERS READ(IN,2001,END=50,ERR=50)STN,SER,TIM,SAT,DAY,HR,XMIN,ERR, $ (ORB2(I),I=1,NVAR),(ORB1(I),I=1,14) IF(ILIST .GE. 1) WRITE(IPR,2003) STN,SER,TIM,SAT,DAY,HR,XMIN,ERR, $ (ORB2(I),I=1,NVAR),(ORB1(I),I=1,14) C DATA CHECKING IF(SER .NE. SERIAL) GO TO 160 IF(TIM .EQ. S) TIM = 1. IF(TIM .EQ. R) TIM = 2. IF(TIM .NE. TCODE) GO TO 165 DLOCK = DAY TLOCK = XMIN + HR * 60. STAT(1) = 1. STAT(2) = 2. ELEM(1) = 14. 105 IF(ORB2(NVAR) .NE. 0.) GO TO 106 NVAR = NVAR - 1 GO TO 105 106 ELEM(2) = NVAR C READ DOPPLER COUNTS K2 = 0 DO 120 I = 1,9 K1 = K2 + 1 K2 = K1 + 3 READ(IN,2002,END=50,ERR=50) (DOP1(K),DOP2(K),K=K1,K2) IF(ILIST .GE. 1) WRITE(IPR,2004) (DOP1(K),DOP2(K),K=K1,K2) C DETECT BLANK CARD = END OF DOPPLER COUNT INPUT DO 110 K = K1,K2 IF(DOP1(K) .NE. 0.) GO TO 120 IF(DOP2(K) .NE. 0.) GO TO 120 110 CONTINUE GO TO 125 120 CONTINUE K2 = K2 + 4 125 NDOP = K2 - 4 DO 126 I = 1,NDOP J = 1 + NDOP - I IF(DOP1(J) .NE. 0. .OR. DOP2(J) .NE. 0.) GO TO 127 126 CONTINUE 127 NDOP = J ELEM(4) = NDOP ELEM(5) = NDOP STAT(4) = 1. STAT(5) = 2. C DIVIDE ALL DOPPLER COUNTS BY TEN (ONE FRACTIONAL DIGIT) DO 130 I = 1,NDOP DOP1(I) = DOP1(I) / 10D0 130 DOP2(I) = DOP2(I) / 10D0 RETURN 150 CONTINUE WRITE(IPR,3001) IFORM STOP C ERROR RETURNS 160 WRITE(IPR,2005) SERIAL,SER STOP 165 WRITE(IPR,2006) TIME,TIM STOP C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE MJVIN ENTERED ) 1001 FORMAT(5I3,5(/7I10)) 1002 FORMAT(8I10) 1003 FORMAT(5X,5I3,5(/5X,7I10)) 1004 FORMAT(5X,12I10) 1005 FORMAT(5X,25HINCONSISTENT RCVR CODE ,2I10) 1006 FORMAT(5X,35HINCONSISTENT NUMBER OF DOPPLERS/MSG ,2I10) 1007 FORMAT(5X,I5,2F12.2) 2001 FORMAT(2X,A2,F4.0,1X,A1,F3.0,F5.0,F4.0,F3.0,1X,10F3.0,4(/1X,7F10.0 $ )) 2002 FORMAT(1X,4(F9.0,F8.0)) 2003 FORMAT(5X,A2,F5.0,1X,A1,4F6.0,5X,10F4.0,4(/10X,7F11.0)) 2004 FORMAT(10X,8F11.0) 2005 FORMAT(5X,35HINCONSISTENT RECEIVER SERIAL NUMBER ,2F10.0) 2006 FORMAT(5X,35HINCONSISTENT RECEIVER TIMING CODE ,2F10.0) 3001 FORMAT(5X,30HWRONG FORMAT CODE IN MJVIN ,I5) END SUBROUTINE SKEP C SCALE KEPLER PARAMETERS IN ORB1 C TP = TIME OF PERIGEE (MIN) UNB00070 C XN = MEAN MOTION (RAD/MIN) UNB00080 C WNOT = ARGUMENT OF PERIGEE (RAD) UNB00090 C WDOT = PRECESSION OF PERIGEE (RAD/MIN) UNB00100 C EXC = ECCENTRICITY UNB00110 C AO = SEMI-MAJOR AXIS (M) UNB00120 C RA = RIGHT ASCENSION OF NODE (RAD) UNB00130 C RADOT = PRECESSION OF NODE (RAD/MIN) UNB00140 C XI = INCLINATION ANGLE (RAD) C GAST = LONGITUDE OF GREENWICH (RAD) UNB00160 C OFFREQ = SATELLITE FREQUENCY OFFSET ( PARTS IN 10(10)) C FG = RECEIVER FREQ OFFSET (PARTS IN 10(10)) C FCA = CLOSEST APPROACH FREQUENCY C SAT = SATELLITE ID NUMBER C DAYINJ = INJECTION DAY C TINJ = INJECTION TIME IN MINUTES C DAYTP = PERIGEE DAY C**MACHINE DEPENDENT STATEMENT UNB00240 IMPLICIT REAL*8(A-H,O-Z) UNB00250 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00260 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION IKEP(14) EQUIVALENCE (FHEAD(55),YEAR),(FHEAD(62),FGO),(FHEAD(63),DFG), $ (FHEAD(64),DAYFGO),(PHEAD(2),SAT),(PHEAD(5),DAYINJ),(PHEAD(6), $ TINJ),(PHEAD(17),STAT),(PHEAD(23),ELEM), $ (ORB1(1),TP),(ORB1(2),XN),(ORB1(3),WNOT), $ (ORB1(4),WDOT),(ORB1(5),EXC),(ORB1(6),AO),(ORB1(7),RA),(ORB1(8), $ RADOT),(ORB1(9),XI),(ORB1(10),GAST),(ORB1(11),DAYTP ),(ORB1(12), $ FCA ),(ORB1(13),FG ),(ORB1(14),OFFREQ) C**MACHINE DEPENDENT STATEMENTS ATAN2(PI,VEL) = DATAN2(PI,VEL) IF(ILIST .GE. 1) WRITE(IPR,1000) UNB00360 IF(STAT .NE. 1. .OR. ELEM .NE. 14.) GO TO 90 I1D7 = 1D7 I1D8 = 1D8 C COPY UNSCALED KEPLERS INTO INTEGER ARRAY FOR DECODING DO 10 I = 1,14 10 IKEP(I) = ORB1(I) C DECODE FIRST DIGIT OF KEPLER PARAMETERS UNB00400 IKEP(1) = MOD(IKEP(1),I1D8) + (IKEP(1)/(4*I1D8))*I1D8 DO 50 I = 2,13 UNB00420 50 IKEP(I) = MOD(IKEP(I),I1D8) * (17-2*(IKEP(I)/I1D8)) IF(ILIST .GE. 3) WRITE(IPR,1002) IKEP C SCALE PARAMETERS UNB00450 TP = IKEP(1) / 1D5 XN = (IKEP(2) / 1D8 + 3D0) WNOT = (IKEP(3) / 1D5) WDOT = (IKEP(4) / 1D8) EXC = IKEP(5) / 1D7 AO = IKEP(6) RA = (IKEP(7) / 1D5) RADOT = (IKEP(8) / 1D8) CI = IKEP(9) / 1D7 SI = IKEP(13) / 1D7 XI = ATAN2(SI,CI) / RADG GAST = (IKEP(10) / 1D5) SAT = MOD(IKEP(11),10000) / 100 DAYINJ = MOD(IKEP(12),10000) / 10 INJT = MOD(IKEP(12),I1D8) / 50000 TINJ = INJT OFFREQ = IKEP(14) / 1D3 C FOR PRE-1973 DATA SAT WAS NOT PART OF MESSAGE - USE AO TO IDENTIFY IF(YEAR .GE. 1973.) GO TO 70 SAT = 0. SAT1 = MOD(IKEP(6),100000) / 10000 IF(SAT1 .EQ. 6.) GO TO 60 IF(SAT1 .EQ. 4.) SAT = 12. IF(SAT1 .EQ. 5.) SAT = 14. 60 SAT1 = MOD(IKEP(6),10000) / 1000 IF(SAT1 .EQ. 2.) SAT = 18. IF(SAT1 .EQ. 3.) SAT = 13. IF(SAT1 .EQ. 4.) SAT = 19. C DETERMINE DAY OF TP UNB00750 70 DAYTP = DAYINJ IF(TP .LT. TINJ) DAYTP = DAYTP + 1. C COMPUTE RECEIVER FREQUENCY OFFSET AND CLOSEST APPROACH FREQ (IN HZ) FG = FGO + DFG * (DAYTP - DAYFGO) FCA = FO * (OFFREQ + FG) * 1D-10 C RESET ORB1 STATUS AND VECTOR LENGTH STAT = 3D0 ELEM = 14D0 IELEM = ELEM IF(ILIST .GE. 3) WRITE(IPR,1001) (ORB1(I),I=1,IELEM) RETURN UNB00780 C ERROR RETURN UNB00790 90 WRITE(IPR,1003) STAT,ELEM STOP C FORMAT STATEMENTS UNB00840 1000 FORMAT(1X,30HSUBROUTINE SKEP ENTERED ) UNB00850 1001 FORMAT(5X,20HSCALED KEPLERS /4(5X,5F20.10/)) 1002 FORMAT(5X,20HUNSCALED KEPLERS ,2(/5X,7I10)) 1003 FORMAT(5X,25HSTATUS OF ORB1 IN SKEP IS ,2F5.1) END SUBROUTINE SEPH C SCALE EPHEMERAL PARAMETERS IN ORB2 C TIME(NVAR) = TIME OF EPHEMERALS (MIN PAST HOUR) UNB00070 C DE(NVAR) = ECCENTRICITY ANOMALY EPHEMERALS (RAD) UNB00080 C DA(NVAR) = SEMI-MAJOR AXIS EPHEMERALS (M) UNB00090 C ETA(NVAR) = OUT OF PLANE EPHEMERALS(M) UNB00100 C IETA = 2 IF EVEN ETA KEPT. = 1 IF ODD ETA KEPT. UNB00110 C**MACHINE DEPENDENT STATEMENT UNB00120 IMPLICIT REAL*8(A-H,O-Z) UNB00130 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00140 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION TIME(17),DE(17),DA(17),ETA(17),VAR(2) DIMENSION IVAR(17) EQUIVALENCE (PHEAD(4),TLOCK),(PHEAD(18),STAT),(PHEAD(24),ELEM), $ (PHEAD(51),XLOCK),(PHEAD(52),ETAI),(ORB2(1),TIME(1)), $ (ORB2(18),DE(1)),(ORB2(35),DA(1)),(ORB2(52),ETA(1)) C**MACHINE DEPENDENT STATEMENT UNB00250 ABS(PI) = DABS(PI) UNB00260 AMOD(PI,VEL) = DMOD(PI,VEL) INT(PI) = IDINT(PI) IF(ILIST .GE. 1) WRITE(IPR,1000) UNB00270 I1D7 = 1D7 I1D8 = 1D8 NVAR = ELEM LOCK = XLOCK IF(STAT .NE. 2.) GO TO 210 IF(NVAR .LT. 9) GO TO 220 IF(NVAR .GT. 17) GO TO 225 DO 5 I = 1,2 5 VAR(I) = 0. X4 = 4D0 X30 = 30D0 X60 = 60D0 X1D1 = 1D1 X1D3 = 1D3 C MOVE ORB2 CONTENTS IF LOCK .LE. 0 IF(LOCK .GT. 0) GO TO 55 ML = 1 - LOCK NVAR = NVAR + ML ELEM = NVAR IF(NVAR .GT. 17) GO TO 225 ML1 = ML + 1 DO 40 I = ML1,NVAR J = ML1 + NVAR - I 40 ORB2(J) = ORB2(J-ML) DO 50 I = 1,ML 50 ORB2(I) = 0. LOCK = 1 UNB00410 XLOCK = LOCK C COPY ORB2 INTO INTEGER ARRAY IVAR FOR DECODING 55 CONTINUE DO 60 I = 1,NVAR 60 IVAR(I) = ORB2(I) IF(ILIST .GE. 3) WRITE(IPR,1001) (IVAR(I),I=1,NVAR) C DECODE FIRST DIGIT AND SEPARATE VARIABLE PARAMETERS UNB00420 DO 100 I = 1,NVAR TIME(I) = (IVAR(I) / (4*I1D8))*10 + MOD(IVAR(I)/I1D7,10) DE(I) = (-1)**(2+IVAR(I)/(2*I1D8))*MOD(IVAR(I)/10000,1000) DA(I) = (-1)**(IVAR(I)/I1D8)*MOD(IVAR(I)/10,1000) ETA(I) = MOD(IVAR(I),10) IF(ILIST .GE. 3 .AND. I .NE. LOCK) WRITE(IPR,1004) TIME(I), $ DE(I),DA(I),ETA(I) IF(ILIST .GE. 3 .AND. I .EQ. LOCK) WRITE(IPR,1005) TIME(I), $ DE(I),DA(I),ETA(I) 100 CONTINUE UNB00500 C FIND A SPAN OF AT LEAST 9 CORRECTLY INCREMENTED VALUES IN TIME VECTOR UNB00510 NCOR = 0 UNB00520 NVAR1 = NVAR - 1 UNB00530 DO 102 I = 1,NVAR1 UNB00540 DT = TIME(I+1) - TIME(I) UNB00550 IF(DT .NE. 1. .AND. DT .NE. -14.) GO TO 101 UNB00560 NCOR = NCOR + 1 UNB00570 IF(NCOR .GE. 8) GO TO 103 UNB00580 GO TO 102 UNB00590 101 NCOR = 0 UNB00600 102 CONTINUE UNB00610 C REJECT PASS IF LT 9 CORRECTLY INCREMENTED TIME VALUES UNB00620 GO TO 220 UNB00630 C REJECT EPHEMERALS IF TIME VALUE IS INCORRECTLY INCREMENTED UNB00640 103 DO 104 J = 1,NVAR UNB00650 DT = TIME(I) + J - I UNB00660 IF(DT .LT. 0.) DT = DT + 15. UNB00670 IF(DT .GT. 14.) DT = DT - 15. UNB00680 IF(TIME(J) .EQ. DT) GO TO 104 UNB00690 TIME(J) = DT UNB00700 DE(J) = 0. UNB00710 DA(J) = 0. UNB00720 ETA(J) = 0. UNB00730 104 CONTINUE UNB00740 C COMBINE OUT OF PLANE DIGITS UNB00750 DO 110 I = 1,NVAR1 UNB00760 SIGN = 1. UNB00770 IF(ETA(I) .LT. 5.) SIGN = -1. UNB00780 ETA1 = 0. UNB00790 IF(ETA(I) .NE. 0.) ETA1 = ABS(ETA(I) - 5.) * 10. UNB00800 110 ETA(I) = SIGN * (ETA1 + ETA(I+1)) UNB00810 ETA(NVAR) = 1D4 C IF LOCKON TIME IS ALREADY KNOWN, COMPARE WITH TIME VECTOR IF(TLOCK .EQ. 0.) GO TO 115 ITLOCK = TLOCK + .1 IF(MOD(ITLOCK,30) .NE. 2*INT(TIME(LOCK)+.1)) GO TO 230 ITIME1 = ITLOCK - 2 * (LOCK - 1) IETA = 1 IF(MOD(ITIME1,4) .NE. 0) IETA = 2 GO TO 135 C COMPUTE VARIANCES OF EVEN AND ODD SETS OF OUT OF PLANE COMPONENTS UNB00830 115 CONTINUE DO 130 I =1,2 UNB00840 N = 0. UNB00850 SUM = 0. UNB00860 SUM2 = 0. UNB00870 DO 120 K = I,NVAR1,2 UNB00880 N = N + 1 UNB00890 SUM = SUM + ETA(K) UNB00900 120 SUM2 = SUM2 + ETA(K) * ETA(K) UNB00910 130 VAR(I) = (SUM2 - SUM**2/N)/(N-1) UNB00920 C DECIDE WHETHER TIME(1) IS IN FIRST OR SECOND HALF OF HOUR UNB00930 IETA = 1 UNB00940 IF(VAR(2) .LT. VAR(1)) IETA = 2 UNB00950 ITIME1 = TIME(1) UNB00960 IF(MOD(IETA + MOD(ITIME1,2),2) .EQ. 0) ITIME1 = ITIME1 + 15 UNB00970 ITIME1 = ITIME1 * 2 UNB00980 135 CONTINUE IF(ILIST .GE. 3) WRITE(IPR,1006) VAR(1),VAR(2),IETA,ITIME1 UNB00990 C SCALE TIME VECTOR,DELETE WRONG VALUES OF ETA, AND REPLACE ALL DELETED C VALUES OF DE,DA,ETA BY 10000 (IMPOSSIBLY LARGE NUMBER) DO 140 I = 1,NVAR UNB01010 TIME(I) = ITIME1 + 2. * (I-1) UNB01020 TIME(I) = AMOD(TIME(I),X60) IF(MOD(IETA + I,2) .EQ. 1) ETA(I) = 1D4 IF(DE(I) .NE. 0. .OR. DA(I) .NE. 0.) GO TO 137 DE(I) = 1D4 DA(I) = 1D4 ETA(I) = 1D4 IF(I .GT. 1) ETA(I-1) = 1D4 137 CONTINUE IF(ILIST .GE. 2) WRITE(IPR,1004) TIME(I),DE(I),DA(I),ETA(I) UNB01050 140 CONTINUE UNB01060 C RESET ORB2 STATUS ETAI = IETA STAT = 4. RETURN UNB01070 C ERROR RETURN UNB01080 210 WRITE(IPR,1002) STAT STOP 220 WRITE(IPR,1008) NVAR NULL = 5 RETURN UNB01120 225 WRITE(IPR,1007) NVAR STOP 230 WRITE(IPR,1009) TLOCK,LOCK,(TIME(I),I=1,NVAR) NULL = 6 RETURN C FORMAT STATEMENTS UNB01130 1000 FORMAT(1X,30HSUBROUTINE SEPH ENTERED ) UNB01140 1001 FORMAT(5X,7I11) 1002 FORMAT(5X,30HWRONG ORB2 STATUS IN SEPH ,F10.1) 1004 FORMAT(10X,4F10.0) UNB01160 1005 FORMAT(1X,4HLOCK,5X,4F10.0) 1006 FORMAT(5X,20H2 VAR + IETA,ITIME1 ,2F10.2,2I10) 1007 FORMAT(5X,25HNVAR TOO LARGE IN SEPH ,I10) 1008 FORMAT(5X,25HNVAR TOO SMALL IN SEPH ,I10) 1009 FORMAT(5X,30HTLOCK AND TIME INCONSISTENT /5X,F10.1,I5/17F4.0) END SUBROUTINE FITEPH C FIT POLYNOMIALS TO ORB2 PARAMETERS, PUTTING COEFFICIENTS IN ORB3 C**MACHINE DEPENDENT STATEMENT UNB00060 IMPLICIT REAL*8(A-H,O-Z) UNB00070 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00080 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION W(1,1) DIMENSION COV(4,4) DIMENSION DE(17),DA(17),ETA(17),F(17),RES(17),BASE(4,17), $ SD(4),CDE(4),CDA(4),CET(4) DIMENSION EPH(17),CEPH(4) EQUIVALENCE(FHEAD(52),POLY),(FHEAD(53),APPSIG),(PHEAD(18),STAT), $ (PHEAD(24),ELEM),(PHEAD(51),XLOCK),(PHEAD(52),ETAI),(ORB1(2),XN) $ ,(ORB2(18),DE(1)),(ORB2(35),DA(1)),(ORB2(52),ETA(1)),(ORB3(1), $ CDE(1)),(ORB3(5),CDA(1)),(ORB3(9),CET(1)) C**MACHINE DEPENDENT STATEMENT UNB00250 SIN(PI) = DSIN(PI) UNB00260 COS(PI) = DCOS(PI) UNB00270 IF(ILIST .GE. 1) WRITE(IPR,1000) UNB00280 IF(STAT .NE. 4.) GO TO 50 LOCK = XLOCK NVAR = ELEM IETA = ETAI NPOLY = POLY DO 40 K = 1,3 DO 10 I = 1,NVAR IF(K .EQ. 1) EPH(I) = DE(I) IF(K .EQ. 2) EPH(I) = DA(I) IF(K .EQ. 3) EPH(I) = ETA(I) 10 CONTINUE C COMPUTE BASE FUNCTIONS J = 0 DO 20 I = 1,NVAR IF(EPH(I) .GT. 1D3) GO TO 20 J = J + 1 F(J) = EPH(I) T = (I - LOCK) * 2D0 BASE(1,J) = 1D0 BASE(2,J) = COS(2. * T * XN * RADG) BASE(3,J) = SIN(2. * T * XN * RADG) BASE(4,J) = T 20 CONTINUE IF(J .LT. NPOLY + 1) GO TO 60 UNB00460 C COMPUTE POLYNOMIAL COEFFICIENTS CEPH CALL LSA(F,J,BASE,NPOLY,4,1,1,W,CEPH,RES,ASIG,COV,ICHK) IF(ICHK .NE. 0) GO TO 80 C CHI-SQUARE TEST OF FIT OF CEPH TO EPH IF(ILIST .GE. 2) WRITE(IPR,1001) (RES(I),I=1,J),ASIG UNB00500 STEST = CHISQ(0.975D0,J-NPOLY) * APPSIG / (J - NPOLY) UNB00510 IF(ASIG .GT. STEST) GO TO 70 UNB00520 DO 30 I = 1,NPOLY IF(K .EQ. 1) CDE(I) = CEPH(I) IF(K .EQ. 2) CDA(I) = CEPH(I) IF(K .EQ. 3) CET(I) = CEPH(I) 30 CONTINUE 40 CONTINUE IF(ILIST .GE. 1) WRITE(6,1004) CDE,CDA,CET UNB00920 RETURN UNB00930 C ERROR RETURNS UNB00940 50 WRITE(IPR,1006) STAT STOP 60 WRITE(IPR,1002) NPOLY,J NULL = 5 RETURN UNB00980 70 WRITE(IPR,1003) STEST,ASIG NULL = 5 RETURN UNB01020 80 WRITE(IPR,1005) NULL = 5 RETURN C FORMAT STATEMENTS UNB01030 1000 FORMAT(1X,30HSUBROUTINE FITEPH ENTERED ) UNB01040 1001 FORMAT(5X,17F7.2) UNB01050 1002 FORMAT(5X,30HNOT ENOUGH VALUES IN FITEPH ,2I10) 1003 FORMAT(5X,30HFAILED CHISQUARE IN FITEPH ,2F15.5) 1004 FORMAT(2X,4F10.4) 1005 FORMAT(5X,25HLSA FAILURE FROM FITEPH ) 1006 FORMAT(5X,30HWRONG ORB2 STATUS IN FITEPH ) END SUBROUTINE TVECT(ITIM) C CONSTRUCT VECTOR OF INTEGRATION TIMES IN ARRAY DOP3 C ITIM = 1 FOR INPUT BR INTEGRATION INTERVALS C 2 FOR INPUT CBR INTEGRATION INTERVALS C 3 FOR INPUT BR INTEGRATION EPOCHS C 4 FOR CBR INPUT INTEGRATION EPOCHS C 5 FOR PROCESSING BR INTEGRATION INTERVALS C 6 FOR PROCESSING CBR INTEGRATION INTERVALS C 7 FOR PROCESSING BR INTEGRATION EPOCHS C 8 FOR PROCESSING CBR INTEGRATION EPOCHS C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), $ DOP7(260),DOP8(260),DOP9(260) DIMENSION TAU1(26),TAU2(26) DIMENSION XS(3) EQUIVALENCE (PHEAD(22),STAT),(PHEAD(26),ELEM1),(PHEAD(28),ELEM2), $ (ORB1(13),FG),(FHEAD(86),TAU1N),(FHEAD(87),TAU1(1)),(FHEAD(113), $ TAU2N),(FHEAD(114),TAU2(1)) IF(ILIST .GE. 1) WRITE(IPR,1000) IF(ITIM .LT. 1 .OR. ITIM .GT. 8) GO TO 20 ZERO = 0D0 C CHOOSE BIT PERIOD BP = BITS IF(MOD(ITIM,2) .EQ. 0) BP = BITR / (1D0 + FG * 1D-10) IF(MOD(ITIM,2) .EQ. 0) CALL SATXYZ(ZERO,XS,S1,3) C CHOOSE INPUT OR PROCESSING TIMES MTAU = TAU1N IF(MOD((ITIM - 1)/4,2) .EQ. 1) MTAU = TAU2N NDOP = ELEM1 T = 0D0 IF(ILIST .GE. 2) WRITE(IPR,1001) NDOP,MTAU DO 10 I = 1,NDOP J = MOD(I-1,MTAU) + 1 TAU = TAU1(J) * BP IF(MOD((ITIM-1)/4,2) .EQ. 1) TAU = TAU2(J) * BP IF(MOD(ITIM,2) .NE. 0) GO TO 5 T1 = T + TAU CALL SATXYZ(T1,XS,S2,3) TAU = TAU - (S2 - S1) / VEL S1 = S2 5 CONTINUE T = T + TAU DOP3(I) = TAU IF(MOD((ITIM-1)/2,2) .EQ. 1) DOP3(I) = T IF(ILIST .GE. 2) WRITE(IPR,1002) I,J,TAU,T,DOP3(I) 10 CONTINUE STAT = 5 + ITIM ELEM2 = ELEM1 RETURN C ERROR RETURN 20 WRITE(IPR,1003) ITIM STOP C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE TVECT ENTERED ) 1001 FORMAT(5X,9HNDOP,NTAU,2I5/14X,1HI,9X,1HJ,7X,3HTAU,12X,1HT, $ 10X,4HDOP3) 1002 FORMAT(5X,2I10,3F15.6) 1003 FORMAT(1X,25HSUBROUTINE TVECT ARGUMENT ,I5) END SUBROUTINE DOPTCA C SCAN 4.6 DOPPLERS TO FIND UNB00040 C TCADOP = TIME OF CLOSEST APPROACH (MIN PAST HOUR) UNB00050 C SLOPE = SLOPE OF DOPPLER CURVE AT CLOSEST APPROACH (M/SEC**2) UNB00060 C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), $ DOP7(260),DOP8(260),DOP9(260) DIMENSION TIME(17) EQUIVALENCE (PHEAD(26),ELEM),(PHEAD(51),XLOCK),(PHEAD(53),TCADOP), $ (PHEAD(54),SLOPE),(PHEAD(55),TCAN),(ORB1(12),FCA),(ORB2(1), $ TIME(1)),(ORB1(13),FG) IF(ILIST .GE. 1) WRITE(IPR,1000) LOCK = XLOCK NDOP = ELEM C FCA = CLOSEST APPROACH FREQUENCY IF(ILIST .GE. 3) WRITE(IPR,1002) FCA C SCAN NONZERO DOPPLERS TO FIND FIRST DOPPLER FREQ HIGHER THAN FCA TIME2 = 0. DO 10 I = 1,NDOP TIME2 = TIME2 + DOP3(I) IF(DOP1(I) .EQ. 0.) GO TO 10 FDOP2 = DOP1(I) / DOP3(I) IF(FDOP2 .LT. FCA) GO TO 10 TIME2 = TIME2 - DOP3(I) / 2. N2 = I IF(ILIST .GE. 3) WRITE(IPR,1003) FDOP2,TIME2,N2 GO TO 20 10 CONTINUE GO TO 80 C SCAN BACKWARDS TO FIND FIRST DOPPLER FREQUENCY LOWER THAN FCA 20 TIME1 = TIME2 + DOP3(N2) / 2. DO 30 J = 1,N2 I = 1 + N2 - J TIME1 = TIME1 - DOP3(I) IF(DOP1(I) .EQ. 0.) GO TO 30 FDOP1 = DOP1(I) / DOP3(I) IF(FDOP1 .GT. FCA) GO TO 30 TIME1 = TIME1 + DOP3(I) / 2. N1 = I IF(ILIST .GE. 3) WRITE(IPR,1004) FDOP1,TIME1,N1 GO TO 40 30 CONTINUE GO TO 80 C COMPUTE SLOPE IN M/SEC/SEC AND TCADOP IN MINUTES UT 40 SLOPE = (FDOP2 - FDOP1) / (TIME2 - TIME1) TCADOP = TIME2 - (FDOP2 - FCA) / SLOPE SLOPE = SLOPE / ((1. + FG*1D-10) * WAVENO) TCADOP = TCADOP / 60. + TIME(LOCK) NTCA = (N1 + N2) / 2 + 0.5 TCAN = NTCA IF(ILIST .GE. 3) WRITE(IPR,1005) SLOPE,TCADOP,NTCA RETURN C ERROR RETURN 80 WRITE(IPR,1001) NULL = 1 RETURN UNB02340 C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE DOPTCA ENTERED ) 1001 FORMAT(5X,20HREJECT-MISSED TCA ) 1002 FORMAT(5X,20HFCA ,F15.3) 1003 FORMAT(5X,20HFDOP2,TIME2,N2 ,F15.3,F15.6,I10) 1004 FORMAT(5X,20HFDOP1,TIME1,N1 ,F15.3,F15.6,I10) 1005 FORMAT(5X,20HSLOPE,TCADOP,NTCA ,F15.3,F15.6,I10) END SUBROUTINE ALERT C COMPARE THEORETICAL AND DOPPLER CURVE TCA,SCA TO DETERMINE CORRECT C ALERT. COMPARE RESULTING LOCKON TIME WITH TLOCK (IF NONZERO). C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION TCA2(10),TEST(10) DIMENSION TIME(17) EQUIVALENCE (PHEAD(3),DLOCK),(PHEAD(4),TLOCK),(PHEAD(5),DAYINJ), $ (PHEAD(6),XINJ),(PHEAD(7),TCA),(PHEAD(8),ECA),(PHEAD(9),ACA), $ (PHEAD(10),DCA),(PHEAD(11),QUAD),(PHEAD(17),STAT),(PHEAD(51), $ XLOCK),(PHEAD(53),TCADOP),(PHEAD(54),SLOPE),(ORB1(11),DAYTP), $ (ORB2(1),TIME(1)) ABS(PI) = DABS(PI) INT(PI) = IDINT(PI) AMOD(PI,VEL) = DMOD(PI,VEL) SIGN(X,Y) = DSIGN(X,Y) IF(ILIST .GE. 1) WRITE(IPR,1000) IF(STAT .NE. 3.) GO TO 110 LOCK = XLOCK X30 = 30D0 X60 = 60D0 XP5 = 0.51D0 C COMPUTE ALL ALERTS FROM TINJ TO TINJ + 800 MINUTES TINJ = XINJ + 1440. * (DAYINJ - DAYTP) K = 0 DO 30 I = 1,10 T = 0. IF(I .EQ. 1) T = TINJ CALL NEXT(T,TCA1,SCA1,ECA1,ACA1,DCA1,QUAD1) IF(TCA1 .GT. TINJ + 800.) GO TO 40 IF(TCA1 .EQ. 0.) GO TO 30 K = K + 1 TCA2(K) = TCA1 ITEST1 = ABS(AMOD(TCA1,X60) - TCADOP) ITEST1 = MIN0(ITEST1,60-ITEST1) ITEST2 = ABS(SCA1 - SLOPE) TEST(K) = ITEST1 + ITEST2 IF(ILIST .GE. 3) WRITE(IPR,1001) K,ITEST1,ITEST2,TEST(K) 30 CONTINUE C FIND WHICH ALERT BEST FITS DOPTCA RESULTS 40 IF(K .EQ. 0) GO TO 90 XMIN = 1D10 DO 50 I = 1,K IF(TEST(I) .GT. XMIN) GO TO 50 XMIN = TEST(I) IFIND = I 50 CONTINUE C REJECT PASS IF NO ALERT FITS WELL IF(XMIN .GT. 5.) GO TO 100 C RECOMPUTE THE CORRECT ALERT TO GET THE GEOMETRICAL PARAMETERS T = TCA2(IFIND) CALL NEXT(T,TCA1,SCA1,ECA1,ACA1,DCA1,QUAD1) TCA = TCADOP+X60*(INT(TCA1/X60)+INT((AMOD(TCA1,X60)-TCADOP)/X30)) ECA = ECA1 ACA = ACA1 DCA = DCA1 QUAD = QUAD1 IF(ILIST .GE. 2) WRITE(IPR,1004) TCA,ECA,ACA,DCA,QUAD C COMPUTE DLOCK1 AND TLOCK1 FROM TCA TLOCK1 = TIME(LOCK) + X60*(INT(TCA/X60) + INT((AMOD(TCA,X60) - $ TIME(LOCK))/X30)) DLOCK1 = DAYTP + INT(TLOCK1/1440.) + INT(SIGN(XP5,TLOCK1) - XP5) TLOCK1 = TLOCK1 - 1440. * (DLOCK1 - DAYTP) C COMPARE WITH INPUT DLOCK AND TLOCK, IF NONZERO IF(DLOCK .EQ. 0.) GO TO 60 IF(ABS(DLOCK-DLOCK1).GT..01.OR.ABS(TLOCK-TLOCK1).GT..01)GO TO 120 RETURN 60 TLOCK = TLOCK1 DLOCK = DLOCK1 IF(ILIST .GE. 3) WRITE(IPR,1007) DLOCK,TLOCK RETURN C ERROR RETURNS 90 WRITE(IPR,1002) NULL = 4 RETURN 100 WRITE(IPR,1003) XMIN NULL = 4 RETURN 110 WRITE(IPR,1005) STAT STOP 120 WRITE(IPR,1006) DLOCK,TLOCK,DLOCK1,TLOCK1 NULL = 6 RETURN C FORMAT STATEMENTS 1000 FORMAT(1X,50HSUBROUTINE ALERT ENTERED ) 1001 FORMAT(5X,15HVISIBLE. TESTS= ,80X,3I5,F5.0) 1002 FORMAT(5X,40HNO VISIBLE PASSES FOR THIS SATELLITE ) 1003 FORMAT(5X,40HBAD BEST FITTING ALERT, WITH TEST VALUE , E20.10) 1004 FORMAT(5X,30HCORRECT TCA,ECA,ACA,DCA,QUAD ,F20.10,4F10.1) 1005 FORMAT(5X,30HWRONG ORB1 STAT IN ALERT , F10.1) 1006 FORMAT(5X,30HINPUT AND ALERT LOCKONS DIFFER ,4F10.1) 1007 FORMAT(5X,20HLOCKON DAY AND TIME ,2F10.1) END SUBROUTINE NEXT(T,TCA,SCA,ECA,ACA,DCA,QUAD) C COMPUTE SATELLITE ALERT NEAREST IN TIME TO T C IF T .EQ. 0. COMPUTE ALERT FOR NEXT PASS AFTER OLD TNODE VALUE C RETURN ZERO OUTPUTS IF NO VISIBLE PASS C TCA = TIME OF CLOSEST APPROACH IN MINUTES C SCA = SLOPE OF DOPPLER CURVE AT CLOSEST APPROACH IN HZ/S C ECA = ELEVATION ANGLE OF SATELLITE AT CLOSEST APPROACH IN DEGREES C ACA = AZIMUTH OF SATELLITE AT CLOSEST APPROACH IN DEGREES C DCA = DIRECTION OF SATELLITE TRAVEL AT CLOSEST APPROACH IN DEGREES C QUAD = QUADRANT OF SATELLITE PASS C 1 = W NS (SATELLITE PASSES TO WEST, TRAVELLING NORTH TO SOUTH) C 2 = W SN C 3 = E NS C 4 = E SN C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION XG(3),X0(3) DIMENSION PLH(3) EQUIVALENCE (FHEAD(14),PLH(1)),(FHEAD(35),RG),(PHEAD(17),STAT), $ (ORB1(1),TP),(ORB1(2),XN),(ORB1(3),WNOT),(ORB1(4),WDOT), $ (ORB1(5),EXC),(ORB1(6),AO),(ORB1(7),RA),(ORB1(8),RADOT),(ORB1(10 $ ),GAST) ABS(PI) = DABS(PI) SIN(PI) = DSIN(PI) COS(PI) = DCOS(PI) ATAN(PI) = DATAN(PI) SQRT(PI) = DSQRT(PI) ATAN2(X,Y) = DATAN2(X,Y) SIGN(X,Y) = DSIGN(X,Y) IF(ILIST .GE. 1) WRITE(IPR,1000) IF(STAT .NE. 3.) GO TO 40 N = 1 IF(T .EQ. 0.) GO TO 5 C INITIALIZE ONLY WHEN T .NE. 0. ONE = 1D0 X60 = 60D0 X90 = 90D0 PI2 = 2D0 * PI RPHI = PLH(1) * RADG XLONG = PLH(2) * RADG SOBS = SIN(RPHI) COBS = COS(RPHI) TOBS = SOBS / COBS REARTH = RG C COMPUTE NODAL CROSSING PARAMETERS PERIOD = 360D0 / (XN - WDOT) WMOTN = (WE - RADOT * RADG) * PERIOD DT = WNOT / (XN + WDOT) TNODE = TP - DT XLNODE = (RA - GAST)*RADG + (WE - RADOT*RADG) * DT IF(XLNODE .LT. 0.) XLNODE = XLNODE + PI2 W2PI = WMOTN / PI2 WPHI = (XN - WDOT) * RADG / X60 WLAM = (WE - RADOT*RADG) / X60 C COMPUTE TNODE,XLNODE IMMEDIATELY PRIOR TO T AN = (T - TNODE) / PERIOD - RPHI / PI2 IF(AN .LT. 0.) N = AN - 0.2 IF(AN .GT. 0.) N = AN + 0.2 C UPDATE TNODE,XLNODE FOR THIS CALL 5 CONTINUE TNODE = TNODE + N * PERIOD XLNODE = XLNODE - N * WMOTN NDAY = XLNODE / PI2 XLNODE = XLNODE - NDAY * PI2 IF(XLNODE .LT. 0.) XLNODE = XLNODE + PI2 IF(ILIST .GE. 3) WRITE(IPR,1001) T,TNODE,XLNODE C INITIALIZE OUTPUT RORBIT = AO COSAMX = REARTH / RORBIT TCA = 0D0 SCA = 0D0 ECA = 0D0 ACA = 0D0 DCA = 0D0 QUAD = 0D0 C COMPUTE X(CA SUBLAT) AND Y (CA SUBLONG) Y0 = XLNODE - XLONG IDIR = SIGN(ONE,COS(Y0)) + 1. X = RPHI IF(IDIR .EQ. 0) X = PI - X DO 10 ITER = 1,10 XX = X Y = Y0 - W2PI * X X = ATAN((TOBS + W2PI * SIN(Y)) / COS(Y)) IF(IDIR .EQ. 0) X = X + PI IF(ABS(X - XX) .LT. 1D-3) GO TO 20 10 CONTINUE C NO CONVERGENCE IF(ILIST .GE. 3) WRITE(IPR,1005) X,Y RETURN C TEST WHETHER PASS IS ABOVE HORIZON 20 SINX = SIN(X) COSX = COS(X) SINY = SIN(Y) COSY = COS(Y) TANY = SINY / COSY CALF = COSX * COSY * COBS + SINX * SOBS IF(CALF .GT. COSAMX) GO TO 30 IF(ILIST .GE. 3) WRITE(IPR,1005) X,Y,CALF RETURN C VISIBLE PASS. COMPUTE TCA,MORE ACCURATE RORBIT,SCA,ECA,ACA,DCA,QUAD 30 TCA = TNODE + PERIOD * X / PI2 RORBIT = AO * (1D0 - EXC * COS(XN*RADG*(TCA-TP))) COSAMX = REARTH / RORBIT SCA = (REARTH / SQRT(1. + (COSAMX**2) - 2.*COSAMX*CALF)) * $ (CALF * (WPHI**2) - COSY *COSX*COBS*(WLAM**2)) SALF = SQRT(1. - CALF**2) ECA = ATAN((CALF - COSAMX)/SALF)/RADG CACA = (SINX - SOBS*CALF)/(COBS*SALF) SACA = SQRT(1. - CACA**2) ACA = ATAN2(SIGN(SACA,TANY),SIGN(CACA,SINX)) / RADG DCA = ACA - SIGN(X90,SINY) ISIDE = SIGN(ONE,TANY) + 1. QUAD = 1 + IDIR / 2 + ISIDE IF(ILIST .GE. 3) WRITE(IPR,1002) TCA,SCA,ECA,ACA,DCA,QUAD RETURN C ERROR RETURN 40 WRITE(IPR,1003) STAT STOP C FORMAT STATEMENTS 1000 FORMAT(1X,50HSUBROUTINE NEXT ENTERED ) 1001 FORMAT(5X,15HT,TNODE,XLNODE ,2F11.5,F10.6) 1002 FORMAT(5X,15HTCA-S-E-A-DCA,Q ,35X,F11.5,F7.3,F5.1,2F7.1,F4.0) 1003 FORMAT(5X,30HWRONG ORB1 STATUS IN NEXT ,F10.1) 1005 FORMAT(5X,15HNONVIS X,Y,CALF ,35X,2F10.6,F8.5) END SUBROUTINE METOBS C OBTAIN SURFACE METEOROLOGICAL OBSERVATIONS UNB00030 C TKO = SURFACE TEMPERATURE IN KELVINS C PSTN = SURFACE PRESSURE IN MB C PWVO = PARTIAL PRESSURE OF WATER VAPOUR IN MB C**MACHINE DEPENDENT STATEMENT UNB00070 IMPLICIT REAL*8(A-H,O-Z) UNB00080 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00090 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) EQUIVALENCE(FHEAD(16),HT),(FHEAD(81),WCODE),(FHEAD(82),TK), $ (FHEAD(83),PSL),(FHEAD(84),PWV),(PHEAD(3),DLOCK),(PHEAD(4), $ TLOCK),(PHEAD(12),TKO),(PHEAD(13),PSTN),(PHEAD(14),PWVO) DATA METDAY/0/,METHR/0/ UNB00150 C**MACHINE DEPENDENT STATEMENT UNB00160 EXP(PI) = DEXP(PI) UNB00170 IF(ILIST .GE. 1) WRITE(IPR,1000) UNB00180 IF(WCODE .EQ. 2.) GO TO 4 C USE STANDARD WEATHER VALUES TKO = TK C = 29.2897 * (TKO + HT / 400.) PSTN = PSL * EXP(-HT / C) PWVO = PWV RETURN C READ OBSERVED WEATHER DATA UNTIL TIME MATCHES TLOCK 4 CONTINUE NBACK = 0 UNB00200 LODAY = DLOCK LOHR = TLOCK / 60D0 5 CONTINUE UNB00210 IF(METDAY .EQ. LODAY .AND. METHR .EQ. LOHR) RETURN UNB00220 TMET = METDAY + METHR / 24D0 UNB00230 TPASS = DLOCK + TLOCK / 1440. IF(TMET .LT. TPASS) GO TO 10 IF(NBACK .GT. 24) GO TO 30 UNB00260 NBACK = NBACK + 1 UNB00270 C BACKSPACE MET C BACKSPACE MET 10 CONTINUE C READ(MET,1001,END=30) METYR,METDAY,METHR,TCO,PSLO,PWVO IF(METHR .LT. 24) GO TO 20 UNB00310 METHR = METHR - 24 UNB00320 METDAY = METDAY + 1 UNB00330 20 CONTINUE UNB00340 TKO = TCO + 273.16 C = 29.2897 * (TKO + HT / 400.) PSTN = PSLO * EXP(-HT / C) GO TO 5 UNB00380 30 CONTINUE UNB00390 WRITE(IPR,1002) METDAY,METHR,LODAY,LOHR UNB00400 RETURN C FORMAT STATEMENTS UNB00460 1000 FORMAT(1X,30HSUBROUTINE METOBS ENTERED ) UNB00470 1001 FORMAT(3I4,3F6.1) UNB00480 1002 FORMAT(10X,24HUNABLE TO MATCH MET TIME ,2I5,12H WITH LOCKON,2I5) UNB00490 END UNB00500 SUBROUTINE POLE C OBTAIN POLE COORDINATES FOR LODAY UNB00320 IMPLICIT REAL*8(A-H,O-Z) BUG***** COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION DPOLE(10),XP(10),YP(10) EQUIVALENCE (PHEAD(5),DAYINJ),(PHEAD(15),XPOLE),(PHEAD(16),YPOLE) $ ,(FHEAD(155),DPOLE(1)),(FHEAD(165),XP(1)),(FHEAD(175),YP(1)) DO 10 I = 1,9 IF(DAYINJ.GE.DPOLE(I) .AND. DAYINJ.LT.DPOLE(I+1)) GO TO 20 10 CONTINUE RETURN 20 XPOLE = XP(I) YPOLE = YP(I) IF(ILIST .GE. 2) WRITE(IPR,1005) XPOLE,YPOLE RETURN 1005 FORMAT(5X,20HPOLE POSNS IN METRES,2F15.5) END SUBROUTINE RANVEL C COMPUTE SATELLITE COORDINATES IN COORDINATE SYSTEM DEFINED BY UNB00030 C SLANT RANGE VECTOR AT CLOSEST APPROACH (X-AXIS) AND VELOCITY UNB00040 C VECTOR AT CLOSEST APPROACH (Y-AXIS) UNB00050 C**MACHINE DEPENDENT STATEMENT UNB00130 IMPLICIT REAL*8(A-H,O-Z) UNB00140 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00150 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION X(3),XDOT(3),R(5),RDOT(3),ROT(3,3),DROT(3,3),NAX(5) DIMENSION CDE(4),CDA(4),CET(4),XG(3) DIMENSION GPROT(3,3) DIMENSION U1(3),U2(3),U3(3) DIMENSION NAX1(3) DATA NAX1(1),NAX1(2),NAX1(3)/3,2,1/ DATA NAX(1),NAX(2),NAX(3),NAX(4),NAX(5)/3,1,3,2,1/ EQUIVALENCE (FHEAD(2),XG(1)),(PHEAD(3),DLOCK),(PHEAD(4),TLOCK), $ (PHEAD(7),TCA),(PHEAD(42),GPROT(1,1)),(ORB1(1),TP),(ORB1(2),XN), $ (ORB1(3),WNOT),(ORB1(4),WDOT),(ORB1(5),EXC),(ORB1(6),AO), $ (ORB1(7),RA),(ORB1(8),RADOT),(ORB1(9),XI),(ORB1(10),GAST), $ (ORB1(11),DAYTP), $ (ORB2(1),CDE(1)),(ORB2(5),CDA(1)),(ORB2(9),CET(1)) $ ,(PHEAD(15),XPOLE),(PHEAD(16),YPOLE) C**MACHINE DEPENDENT STATEMENT UNB00290 ARCOS(PI) = DARCOS(PI) UNB00300 ATAN2(PI,VEL) = DATAN2(PI,VEL) UNB00310 COS(PI) = DCOS(PI) UNB00320 SIN(PI) = DSIN(PI) UNB00330 SQRT(PI) = DSQRT(PI) UNB00340 IF(ILIST .GE. 1) WRITE(IPR,1000) UNB00350 T = TCA - TLOCK + 1440. * (DAYTP - DLOCK) C COMPUTE UNIT SLANT RANGE VECTOR U1 CT = COS(2. * T * XN * RADG) ST = SIN(2. * T * XN * RADG) DE = (CDE(1) + CDE(2) * CT + CDE(3) * ST + CDE(4) * T) * RADG/1D4 DA = (CDA(1) + CDA(2) * CT + CDA(3) * ST + CDA(4) * T) * 1D1 ET = (CET(1) + CET(2) * CT + CET(3) * ST + CET(4) * T) * 1D1 DT = TCA - TP XE = XN * RADG * DT + EXC * SIN(XN * RADG * DT) + DE XA = AO + DA X(1) = XA * ( COS(XE) - EXC) UNB00590 X(2) = XA * SIN(XE) UNB00600 X(3) = ET R(1) = -(WNOT - WDOT * DT) * RADG R(2) = - XI * RADG R(3) = -(RA + RADOT * DT) * RADG + GAST * RADG + WE * DT R(4) = - XPOLE / 6356D3 R(5) = -YPOLE / 6356D3 CALL ROTREF(5,NAX,R,ROT) DO 20 I = 1,3 U1(I) = 0. DO 10 J = 1,3 10 U1(I) = U1(I) + ROT(I,J) * X(J) 20 U1(I) = U1(I) - XG(I) RU1 = SQRT(U1(1)**2 + U1(2)**2 + U1(3)**2) DO 30 I = 1,3 30 U1(I) = U1(I) / RU1 C COMPUTE UNIT VELOCITY VECTOR U2 DEDOT = (2.*XN*RADG*(-CDE(2)*ST + CDE(3)*CT) + CDE(4)) * RADG/1D4 DADOT = (2.*XN*RADG*(-CDA(2)*ST + CDA(3)*CT) + CDA(4)) * 1D1 ETDOT = (2.*XN*RADG*(-CET(2)*ST + CET(3)*CT) + CET(4)) * 1D1 XEDOT = XN*RADG*(1D0 + EXC * COS(XN*RADG*DT)) + DEDOT XDOT(1) = DADOT * ( COS(XE) - EXC) - XEDOT * XA * SIN(XE) XDOT(2) = DADOT * SIN(XE) + XEDOT * XA * COS(XE) XDOT(3) = ETDOT UNB00670 RDOT(1) = WDOT * RADG RDOT(2) = 0. UNB00770 RDOT(3) = -(RADOT*RADG - WE) CALL ROTDOT(3,NAX,R,RDOT,DROT) UNB00790 DO 40 I = 1,3 U2(I) = 0. DO 40 J = 1,3 40 U2(I) = U2(I) + ROT(I,J) * XDOT(J) + DROT(I,J) * X(J) RU2 = SQRT(U2(1)**2 + U2(2)**2 + U2(3)**2) DO 50 I = 1,3 50 U2(I) = U2(I) / RU2 C COMPUTE UNIT VECTOR U3 = U1 X U2 DO 60 I = 1,3 J = MOD(I,3) + 1 K = MOD(J,3) + 1 60 U3(I) = U1(J) * U2(K) - U1(K) * U2(J) C U1 U2 AND U3 ARE THE ROW VECTORS OF TRANSFORMATION MATRIX GPROT DO 70 I = 1,3 GPROT(1,I) = U1(I) GPROT(2,I) = U2(I) 70 GPROT(3,I) = U3(I) IF(ILIST .GE. 3) WRITE(IPR,1003) GPROT C ALTERNATIVE ALGORITHM R(1) = ATAN2(U1(2),U1(1)) R(2) = -ATAN2(U1(3),SQRT(U1(1)**2+U1(2)**2)) CALL ROTREF(2,NAX1,R,ROT) DO 80 I = 1,3 U3(I) = 0. DO 80 J = 1,3 80 U3(I) = U3(I) + ROT(I,J) * U2(J) R(3) = ATAN2(U3(3),U3(2)) ALF = ATAN2(SQRT(U3(2)**2 + U3(3)**2),U3(1))/RADG CALL ROTREF(3,NAX1,R,GPROT) IF(ILIST .GE. 3) WRITE(IPR,1003) GPROT,ALF RETURN UNB01180 C FORMAT STATEMENTS UNB01190 1000 FORMAT(1X,30HSUBROUTINE RANVEL ENTERED ) UNB01200 1003 FORMAT(5X,3E20.10) END SUBROUTINE SATXYZ(SEC,XS,S,ICS) C COMPUTE SATELLITE COORDINATES XS FOR SEC SECONDS AFTER LOCKON C ICS = COORDINATE SYSTEM CODE FOR SYSTEM TO WHICH XS REFERRED C 1 ORBITAL COORDINATE SYSTEM C 2 AVERAGE TERRESTRIAL SYSTEM C 3 TOPOCENTRIC COORDINATE SYSTEM C 4 GUIER PLANE COORDINATE SYSTEM C S = SATELLITE RADIUS VECTOR IF ICS = 1,2 C SLANT RANGE IF ICS = 3,4 C**MACHINE DEPENDENT STATEMENT UNB00070 IMPLICIT REAL*8(A-H,O-Z) UNB00080 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00090 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION R(5),ROT(3,3),NAX(5),NAX1(3),X(3),XS(3),CDE(4),CDA(4), $ CET(4),XGP(3) DIMENSION XG(3),TOPROT(3,3),GPROT(3,3) EQUIVALENCE (ORB1(1),TP),(ORB1(2),XN),(ORB1(3),WNOT),(ORB1(4), $ WDOT),(ORB1(5),EXC),(ORB1(6),AO),(ORB1(7),RA),(ORB1(8),RADOT), $ (ORB1(9),XI),(ORB1(10),GAST),(ORB1(11),DAYTP),(FHEAD(2),XG(1)), $ (FHEAD(26),TOPROT(1,1)),(PHEAD(3),DLOCK),(PHEAD(4),TLOCK), $ (PHEAD(15),XPOLE),(PHEAD(16),YPOLE),(PHEAD(17),STAT),(PHEAD(42), $ GPROT(1,1)),(ORB3(1),CDE(1)),(ORB3(5),CDA(1)),(ORB3(9),CET(1)) $ ,(PHEAD(30),XGP(1)) DATA NAX(1),NAX(2),NAX(3),NAX(4),NAX(5)/3,1,3,2,1/ UNB00260 C**MACHINE DEPENDENT STATEMENT UNB00270 SIN(PI) = DSIN(PI) UNB00280 COS(PI) = DCOS(PI) UNB00290 SQRT(PI) = DSQRT(PI) IF(ILIST .GE. 1 .AND. SEC .EQ. 0.) WRITE(IPR,1000) C COMPUTE EPHEMERAL PARAMETERS AT T T = SEC / 60D0 CT = COS(2. * T * XN * RADG) ST = SIN(2. * T * XN * RADG) DE = (CDE(1) + CDE(2) * CT + CDE(3) * ST + CDE(4) * T) * RADG/1D4 DA = (CDA(1) + CDA(2) * CT + CDA(3) * ST + CDA(4) * T) * 1D1 ET = (CET(1) + CET(2) * CT + CET(3) * ST + CET(4) * T) * 1D1 C COMPUTE SATELLITE COORDINATES IN ORBITAL COORDINATE SYSTEM UT = T + TLOCK DT = UT - TP + (DLOCK - DAYTP) * 1440. XE = XN * RADG * DT + EXC * SIN(XN * RADG * DT) + DE XA = AO + DA XS(1) = XA * (COS(XE) - EXC) XS(2) = XA * SIN(XE) XS(3) = ET S = SQRT(XS(1)**2 + XS(2)**2 + XS(3)**2) IF(ICS .GT. 1) GO TO 10 IF(ILIST .GE. 4) WRITE(IPR,1001) SEC,XS,S RETURN C COMPUTE ROTATION MATRIX FROM ORBITAL TO AVERAGE TERRESTRIAL SYSTEM 10 CONTINUE R(1) = -(WNOT - WDOT * DT) * RADG R(2) = - XI * RADG R(3) = -(RA + RADOT * DT) * RADG + GAST * RADG + WE * DT R(4) = - XPOLE / 6356D3 R(5) = -YPOLE / 6356D3 CALL ROTREF(5,NAX,R,ROT) C ROTATE ORBITAL COORDINATES INTO AVERAGE TERRESTRIAL COORDINATES DO 20 K = 1,3 X(K) = XS(K) 20 XS(K) = 0D0 DO 30 K = 1,3 DO 30 L = 1,3 30 XS(K) = XS(K) + ROT(K,L) * X(L) IF(ICS .GT. 2) GO TO 40 IF(ILIST .GE. 4) WRITE(IPR,1002) SEC,XS,S RETURN C TRANSLATE FROM GEOCENTRE TO TOPOCENTRE 40 DO 50 K = 1,3 X(K) = XS(K) - XG(K) 50 XS(K) = 0D0 S = SQRT(X(1)**2 + X(2)**2 + X(3)**2) ICS1 = ICS - 2 GO TO (60,80),ICS1 C ROTATE TRANSLATED AVERAGE TERRESTRIAL COORDINATES INTO TOPOCENTRIC SYST 60 DO 70 K = 1,3 DO 70 L = 1,3 70 XS(K) = XS(K) + TOPROT(K,L) * X(L) IF(ILIST .GE. 4) WRITE(IPR,1003) SEC,XS,S RETURN C ROTATE TRANSLATED AVERAGE TERRESTRIAL COORDINATES INTO GUIER PLANE 80 DO 90 K = 1,3 DO 90 L = 1,3 90 XS(K) = XS(K) + GPROT(K,L) * X(L) S = SQRT((XS(1)-XGP(1))**2+(XS(2)-XGP(2))**2+(XS(3)-XGP(3))**2) IF(ILIST .GE. 4) WRITE(IPR,1004) SEC,XS,S RETURN C FORMAT STATEMENTS UNB00760 1000 FORMAT(1X,30HSUBROUTINE SATXYZ ENTERED ) UNB00770 1001 FORMAT(5X,25HORBITAL COORDINATES ,F15.6,4F15.3) 1002 FORMAT(5X,25HAVG TERREST COORDINATES ,F15.6,4F15.3) 1003 FORMAT(5X,25HTOPOCENTRIC COORDINATES ,F15.6,4F15.3) 1004 FORMAT(5X,25HGUIER PLANE COORDINATES ,F15.6,4F15.3) END SUBROUTINE COMBIN C COMBINE DATA VECTORS IN ARRAY DUMMY FROM INPUT INTERVALS TO PROCESSING C INTERVALS C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) DIMENSION TAUN(15) DIMENSION STAT(3),ELEM(3),DOPA(7) EQUIVALENCE (PHEAD(20),STAT(1)),(PHEAD(26),ELEM(1)) $ ,(FHEAD(86),TAU1N),(FHEAD(113),TAU2N),(FHEAD(140),TAUN(1)) $ ,(PHEAD(55),TCAN) IF(ILIST .GE. 1) WRITE(IPR,1000) IF(STAT(1) .NE. 1. .OR. STAT(2) .NE. 2.) GO TO 65 IF(ELEM(1) .NE. ELEM(2)) GO TO 70 NTCA = TCAN NDOP = ELEM(1) NTAU1 = TAU1N NTAU2 = TAU2N NV = 2 IF(ILIST .GE. 2) WRITE(IPR,1001) NTAU1,NTAU2,(TAUN(I),I=1,NTAU2) IF(ILIST .GE. 3) WRITE(IPR,1006) C N = INPUT INTERVAL POINTER (.LE. NDOP). M = PROCESSING INTER POINTER C COMPUTE NTAU2 PROCESSING INTERVALS/TWO MINUTE MESSAGE, LOOPING OVER C MAXIMUM OF 10 MESSAGES UNTIL NDOP REACHED N = 0 M = 0 DO 30 J = 1,10 DO 30 K = 1,NTAU2 M = M + 1 DO 10 I = 1,NV 10 DOPA(I) = 0D0 C COMBINE L2 INPUT INTERVALS INTO CURRENT PROCESSING INTERVAL L2 = TAUN(K) IZERO = 0 DO 20 L = 1,L2 N = N + 1 IF(ILIST.GE.3) WRITE(IPR,1002)N,(DUMMY(N,I),I=1,NV) IF(N .GT. NDOP) GO TO 40 IF(N .EQ. NTCA) TCAN = M IF(DUMMY(N,1).EQ.0. .OR. DUMMY(N,2).EQ.0.) IZERO = 1 DO 15 I = 1,NV 15 DOPA(I) = DOPA(I) + DUMMY(N,I) 20 CONTINUE 25 DO 27 I = 1,NV IF(IZERO .NE. 0) DOPA(I) = 0D0 27 DUMMY(M,I) = DOPA(I) IF(ILIST.GE.3) WRITE(IPR,1005) M,(DUMMY(M,I),I=1,NV) 30 CONTINUE 40 M = M - 1 NDOP = M DO 60 I = 1,3 60 ELEM(I) = NDOP RETURN C ERROR RETURNS 65 WRITE(IPR,1003) STAT STOP 70 WRITE(IPR,1004) ELEM STOP C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE COMBIN ENTERED ) 1001 FORMAT(5X,50HNUMBER OF INPUT AND PROCESS INTERVALS/TWO MINUTES , $ 2I10/5X,15F4.0) 1002 FORMAT(5X,I5,6F13.2) 1003 FORMAT(5X,30HWRONG DOP STATUS IN COMBIN ,3F10.1) 1004 FORMAT(5X,30HWRONG DOP ELEM IN COMBIN ,3F10.1) 1005 FORMAT(40X,I5,6F13.2) 1006 FORMAT(15X,10HINPUT DATA,20X,15HPROCESSING DATA/16X,4HDOP1,6X, $ 4HDOP2,16X,4HDOP1,6X,4HDOP2) END SUBROUTINE EDIT C INPUT C DOP1 = HIGH CHANNEL DOPPLERS C DOP2 = LOW CHANNEL DOPPLERS C OUTPUT C DOP7 = ELEVATION ANGLES C DOP8 = EDIT STRING C EDIT CODES C 1 = ACCEPTED DOPPLER C 0 = DOPPLER NOT OBSERVED C -1 = LOSS OF LOCK C -2 = HIGH CHANNEL OUTSIDE LIMITS C -3 = LOW CHANNEL OUTSIDE LIMITS C -4 = ELEVATION ANGLE TOO LOW C -5 = DOPPLER CURVE SLOPE TOO SMALL C -6 = DOPPLER INTERVAL TOO FAR FROM CLOSEST APPROACH C -7 = SYMMETRY EDIT C -8 = COMMON EDIT (MULTISTATION SOLUTION) C -9 = RESIDUAL TOO HIGH C -10 = MISCLOSURE TOO HIGH C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), $ DOP7(260),DOP8(260),DOP9(260) DIMENSION XS(3) EQUIVALENCE (FHEAD(67),ELCUT),(FHEAD(68),CUT),(FHEAD(69),SLOPEM), $ ,(FHEAD(73),DOPMIN) BUG***** $ (PHEAD(26),ELEMD),(PHEAD(55),TCAN) BUG***** ATAN(PI) = DATAN(PI) SQRT(PI) = DSQRT(PI) IF(ILIST .GE. 1) WRITE(IPR,1000) MINDOP = DOPMIN NOSYM = 1 NOSYM = 0 NDOP = ELEMD NDOP1 = NDOP + 1 NTCA = TCAN NCUT = CUT NSYM = MIN0(NTCA-1,NDOP-NTCA) IF(ILIST .GE. 2) WRITE(IPR,1003) NDOP,NDOP1,NTCA,NCUT,NSYM IF(ILIST .GE. 3) WRITE(IPR,1004) C FILL DOP7 WITH ELEVATION ANGLES T = 0. CALL SATXYZ(T,XS,S,3) DOP7(1) = ATAN(XS(3)/SQRT(XS(1)**2 + XS(2)**2)) / RADG DO 10 I = 1,NDOP T = T + DOP3(I) CALL SATXYZ(T,XS,S,3) 10 DOP7(I+1) = ATAN(XS(3)/SQRT(XS(1)**2 + XS(2)**2)) / RADG C FILL DOP8 WITH EDITING CODES DO 20 I = 1,NDOP1 IF(DOP1(I) .GT. 0. .OR. DOP2(I) .GT. 0.) DOP8(I) = 1. IF((DOP7(I) .LT. ELCUT .OR. DOP7(I+1) .LT. ELCUT) .AND. $ DOP8(I) .GT. 0.) DOP8(I) = -4. IF(I .LT. 2 .OR. I .GT. NDOP) GO TO 15 IF((DOP1(I)/DOP3(I)-DOP1(I-1)/DOP3(I-1))*2./(DOP3(I)+DOP3(I-1)) $ .LT. SLOPEM .AND. DOP8(I) .GT. 0.) DOP8(I) = -5. 15 CONTINUE IF(IABS(I-NTCA) .GT. NCUT .AND. DOP8(I) .GT. 0.) DOP8(I) = -6. 20 CONTINUE IF(NOSYM .EQ. 1) GO TO 35 DO 30 I = 1,NSYM IF(DOP8(NTCA+I) .GT. 0. .AND. DOP8(NTCA-I) .GT. 0.) GO TO 30 IF(DOP8(NTCA-I) .GT. 0.) DOP8(NTCA-I) = -7. IF(DOP8(NTCA+I) .GT. 0.) DOP8(NTCA+I) = -7. 30 CONTINUE 35 CONTINUE IF(ILIST .GE. 3) WRITE(IPR,1001) (I,DOP7(I),DOP8(I),I=1,NDOP1) N = 0 DO 40 I = 1,NDOP IF(DOP8(I) .GT. 0.) N = N + 1 40 CONTINUE IF(ILIST .GE. 1) WRITE(IPR,1002) N IF(N .LT. MINDOP) GO TO 50 RETURN C ERROR RETURNS 50 WRITE(IPR,1005) N NULL = 2 RETURN 1000 FORMAT(1X,30HSUBROUTINE EDIT ENTERED ) 1001 FORMAT(5X,I5,F5.1,F5.0) 1002 FORMAT(5X,20HNUMBER OF DOPS ,I5) 1003 FORMAT(5X,6I5) 1004 FORMAT(9X,11HI ELEV EDIT ) 1005 FORMAT(5X,11HREJECT-ONLY ,I5,5H DOPS ) END SUBROUTINE IONOS C IONOSPHERIC REFRACTION CORRECTION OF DOPPLER DATA C INPUT C DOP1 = HIGH CHANNEL DOPPLERS C DOP2 = LOW CHANNEL DOPPLERS C DOP3 = INTEGRATION INTERVALS C OUTPUT C DOP4 = VACUUM DOPPLERS C DOP5 = IONOSPHERIC CORRECTION C DOP8 = EDIT STRING C**MACHINE DEPENDENT STATEMENT UNB00090 IMPLICIT REAL*8(A-H,O-Z) UNB00100 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00110 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), $ DOP7(260),DOP8(260),DOP9(260) EQUIVALENCE(PHEAD(20),STAT1),(PHEAD(21),STAT2),(PHEAD(22),STAT3), $ (PHEAD(26),ELEM1),(PHEAD(27),ELEM2),(PHEAD(28),ELEM3), $ (ORB1(12),FCA) C**MACHINE DEPENDENT STATEMENT UNB00210 ABS(PI) = DABS(PI) UNB00220 IF(ILIST .GE. 1) WRITE(IPR,1000) UNB00230 IF(STAT1 .NE. 1. .OR. STAT2 .NE. 2.) GO TO 90 IF(ELEM1 .NE. ELEM2 .OR. ELEM1 .NE. ELEM3) GO TO 100 NDOP = ELEM1 IF(ILIST .GE. 3) WRITE(IPR,1004) C FMAX = MAXIMUM DOPPLER SHIFT FREQUENCY ( IN HZ) C REFMAX = MAXIMUM REFRACTION FREQUENCY (IN HZ) FMAX = 8000. REFMAX = 10. DO 20 I = 1,NDOP DOP4(I) = 0D0 FDOP = 0D0 FREF = 0D0 DOP5(I) = REFMAX * DOP3(I) C CHECK FOR LOSS OF LOCK IF((DOP1(I).EQ.0. .OR. DOP2(I).EQ.0.).AND.DOP8(I).GT.0.)DOP8(I) $ = -1. C CHECK IF HIGH CHANNEL OUTSIDE LIMITS IF((ABS(DOP1(I)/DOP3(I)-FCA).GT.FMAX).AND.DOP8(I).GT.0.) $ DOP8(I) = -2. C CHECK IF LOW CHANNEL OUTSIDE LIMITS AND CONVERT TO REFRACTION COUNT C MARCONI/MAGNAVOX FORMAT IF(ABS(DOP1(I) - DOP2(I))*(9D0/55D0)/DOP3(I) .LT. REFMAX) $ DOP5(I) = (DOP1(I) - DOP2(I)) * 9D0 / 55D0 C ITT FORMAT IF(ABS(DOP2(I) - 2000.)*(24D0/55D0)/DOP3(I) .LT. REFMAX) $ DOP5(I) = -(DOP2(I) - 2000D0) * 24D0 / 55D0 IF(ABS(DOP5(I)).GE.REFMAX*DOP3(I).AND.DOP8(I).GT.0.) DOP8(I)=-3. IF(DOP8(I) .EQ. -3.) DOP5(I) = 0D0 IF(DOP1(I) .GT. 0.) DOP4(I) = DOP1(I) + DOP5(I) IF(DOP4(I) .GT. 0.) FDOP = DOP4(I) / DOP3(I) - FCA IF(DOP8(I) .GT. 0.) FREF = DOP5(I) / DOP3(I) IF(ILIST .GE. 3) WRITE(IPR,1001) I,DOP1(I),DOP2(I),DOP3(I), $ DOP4(I),DOP5(I),FDOP,FREF,DOP8(I) 20 CONTINUE RETURN C ERROR RETURNS 90 WRITE(IPR,1002) STAT1,STAT2,STAT3 STOP 100 WRITE(IPR,1003) ELEM1,ELEM2,ELEM3 STOP C FORMAT STATEMENTS UNB00590 1000 FORMAT(1X,30HSUBROUTINE IONOS ENTERED ) UNB00600 1001 FORMAT(5X,I5,2F13.2,F13.6,4F13.3,F5.0) 1002 FORMAT(5X,25HWRONG DOP STATUS IN IONOS , 3F10.1) 1003 FORMAT(5X,35HDOP ARRAY LENGTHS DIFFER IN IONOS ,3F10.1) 1004 FORMAT(9X,1HI,6X,4HDOP1,9X,4HDOP2,9X,4HDOP3,9X,4HDOP4,9X,4HDOP5, $ 9X,4HFDOP,9X,9HFREF DOP8 ) END UNB00640 SUBROUTINE TROPO C COMPUTE HOPFIELD TWO-QUARTIC TROPOSPHERIC REFRACTION CORRECTION UNB00070 C COMPUTE TROPOSPHERIC REFRACTION CORRECTION TO SLANT RANGE DRHO UNB01040 C INPUT PARAMETERS UNB01050 C TK = SURFACE TEMPERATURE IN KELVINS UNB01060 C P = TOTAL SURFACE PRESSURE IN MILLIBARS UNB01070 C PWV = PARTIAL PRESSURE OF WATER VAPOUR IN MILLIBARS UNB01080 C DOP4 = DOPPLER COUNTS C DOP7 = ELEVATION ANGLES C OUTPUT C DOP4 = CORRECTED DOPPLER COUNTS C DOP6 = TROPOSPHERIC CORRECTIONS C**MACHINE DEPENDENT STATEMENT UNB00080 IMPLICIT REAL*8(A-H,O-Z) UNB00090 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00100 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), $ DOP7(260),DOP8(260),DOP9(260) COMMON /CTROP/ XNT(2),H0(2) EQUIVALENCE (PHEAD(12),TKO),(PHEAD(13),PSTN),(PHEAD(14),PWVO) $ ,(PHEAD(26),DOPNUM) IF(ILIST .GE. 1) WRITE(IPR,1000) UNB00390 NDOP = DOPNUM XNT(1) = 77.6 * PSTN / TKO XNT(2) = 77.6 * 4810D0 * PWVO / TKO**2 H0(1) = 40136D0 + 148.72D0 * (TKO - 273.16) H0(2) = 11000D0 IF(ILIST .GE. 2) WRITE(IPR,1005) TKO,PSTN,PWVO,XNT,H0 IF(ILIST .GE. 2) WRITE(IPR,1003) ELEV = DOP7(1) DRHO1 = YION(ELEV) DO 100 I = 1,NDOP ELEV = DOP7(I+1) DRHO2 = YION(ELEV) DOPOBS = DOP4(I) IF(ELEV .GT. 0.) DOP6(I) = - WAVENO * (DRHO2 - DRHO1) IF(DOP4(I) .GT. 0.) DOP4(I) = DOP4(I) + DOP6(I) IF(ILIST .GE. 2) WRITE(IPR,1002) I,DOPOBS,DOP6(I),DOP4(I), $ DRHO1,DRHO2 DRHO1 = DRHO2 100 CONTINUE RETURN UNB01910 C FORMAT STATEMENTS UNB02010 1000 FORMAT(1X,30HSUBROUTINE TROPO ENTERED ) UNB02020 1002 FORMAT(5X,I5,5F13.3) 1005 FORMAT(5X,20HTEMP,PRESS,WVPRESS ,3F10.2/5X, $ 20HSURFACE REFRACTIVITY ,2F10.2/5X,20HSCALE HEIGHTS , $ 2F10.2) 1003 FORMAT(9X,1HI,9X,30HDOP4 + DOP6 = DOP4 ,8X,5HDRHO1, $ 8X,5HDRHO2) END UNB02090 FUNCTION YION(ELEV) C YIONOULIS ALGORITHM FOR HOPFIELD TROPOSPHERIC RANGE CORRECTION C HT = OBSERVERS ORTHOMETRIC HEIGHT IN M UNB01090 C RG = OBSERVERS RADIAL DISTANCE FROM GEOCENTRE UNB01100 C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DUMMY(260,9) COMMON /CTROP/ XNT(2),H0(2) DIMENSION DR(2),EL(2),RTRO(2),HTRO(2),G(4),C3(4),NLO(2) DATA EL(1),EL(2),NLO(1),NLO(2)/17D0,7D0,2,1/ DATA C3(1),C3(2),C3(3),C3(4),G(1),G(2),G(3),G(4)/1D0,3D0,3D0,1D0, UNB00310 * 1D0,.25D0,.125D0,.078125D0/ UNB00320 EQUIVALENCE(FHEAD(16),HT),(FHEAD(35),RG) C**MACHINE DEPENDENT STATEMENT SQRT(PI) = DSQRT(PI) COS(PI) = DCOS(PI) SIN(PI) = DSIN(PI) C INITIALIZE HTRO(1) = H0(1) - HT UNB01170 HTRO(2) = H0(2) - HT UNB01180 RTRO(1) = RG + HTRO(1) UNB01190 RTRO(2) = RG + HTRO(2) UNB01200 XL1 = RG * SIN(ELEV * RADG) XL2 = RG * COS(ELEV * RADG) DO 80 K = 1,2 UNB01370 W1 = RTRO(K) + XL2 UNB01380 W2 = RTRO(K) - XL2 UNB01390 RW = SQRT(W1*W2) UNB01400 W2W1 = W2 / W1 UNB01410 W1W2M1 = W1 / W2 - 1. UNB01420 HW2 = HTRO(K) / W2 UNB01430 IF(ELEV .LT. EL(K)) GO TO 50 C HIGH ELEVATION ANGLE ALGORITHM UNB01450 SUMP = 0. UNB01460 DO 40 IP = 1,3 UNB01470 SUMN = 0. UNB01480 DO 30 N = 1,IP UNB01490 30 SUMN = SUMN + G(N)*G(IP-N+1)*W2W1**N UNB01500 40 SUMP = SUMP + (1./(IP+5.))*HW2**(IP+1)*(2.*G(IP+1)*(1.+ UNB01510 * W2W1**(IP+1))-SUMN) UNB01520 DR(K) = XNT(K)*1D-6*(RW-XL1-0.8*HTRO(K)*RTRO(K)/RW-RW*SUMP) UNB01530 GO TO 80 UNB01540 C LOW ELEVATION ANGLE ALGORITHM UNB01550 50 SUMN = 0. UNB01560 NP = NLO(K) UNB01570 DO 70 N = 1,4 UNB01580 SUMP = 0. UNB01590 DO 60 IP = 1,NP UNB01600 C1 = 1. + 2. * (IP + N) UNB01610 60 SUMP = SUMP+(-1.)**(IP-1)*G(IP)/C1*W1W2M1**(-IP)*(1.-(1.UNB01620 * -HW2)**(C1/2.)) UNB01630 C2 = (2. * N + 1.) / 2. UNB01640 70 SUMN = SUMN+(-1.)**(N-1)*C3(N)*(1./C2*(1.-(1.-HW2)**C2)+ UNB01650 * SUMP) UNB01660 DR(K) = XNT(K)*1D-6*(-XL1+4.*W2*HW2**(-4)* SQRT(W1W2M1)*SUMNUNB01670 * ) UNB01680 80 CONTINUE UNB01690 YION = DR(1) + DR(2) RETURN END SUBROUTINE GUIER C COMPUTE TWO DIMENSIONAL FIX ON GUIER PLANE (PLANE DEFINED BY SLANT UNB00030 C RANGE AND VELOCITY VECTORS AT CLOSEST APPROACH) UNB00040 C**MACHINE DEPENDENT STATEMENT UNB00090 IMPLICIT REAL*8(A-H,O-Z) UNB00100 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST UNB00110 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), $ DOP7(260),DOP8(260),DOP9(260) DIMENSION DX(3),BASE(3,260),F(260),C(3),R(260),W(1,1),X(3), $ XS(3),COV(3,3),SD(3) DIMENSION XAT(3),GPROT(3,3),COVGP(3,3),COVAT(3,3) DIMENSION XPLH(3),TOPROT(3,3),COVPL(3,3) DIMENSION SXYZ(3,3),WT(3,3) DIMENSION KWAD(4) DATA KWAD(1),KWAD(2),KWAD(3),KWAD(4)/4HW NS,4HW SN,4HE NS,4HE SN/ EQUIVALENCE (FHEAD(58),DELAY),(FHEAD(60),TCODE),(FHEAD(61),DOPSIG) $ ,(FHEAD(70),YMISCL),(FHEAD(71),ALF1),(FHEAD(73),DOPMIN), $ (FHEAD(74),XITER),(FHEAD(75),EPS1),(FHEAD(76),EPS2),(FHEAD(77), $ EPS3),(FHEAD(78),EPS4),(FHEAD(79),EPS5),(FHEAD(80),ALF2), $ (PHEAD(2),SAT),(PHEAD(3),DLOCK),(PHEAD(4),TLOCK),(PHEAD(5), $ DAYINJ),(PHEAD(6),TINJ),(PHEAD(8),ECA),(PHEAD(10),DCA), $ (PHEAD(11),QUAD),(PHEAD(26),ELEMD),(PHEAD(29),COCODE),(PHEAD(30) $ , X(1)),(PHEAD(33),COV(1,1)),(ORB1(14),OFFREQ),(ORB1(13),FG) $ ,(PHEAD(42),GPROT(1,1)) $ ,(FHEAD(26),TOPROT(1,1)) $ ,(FHEAD(5),SXYZ(1,1)),(FHEAD(65),FGOSIG) C**MACHINE DEPENDENT STATEMENT UNB00300 ABS(PI) = DABS(PI) UNB00310 SQRT(PI) = DSQRT(PI) UNB00320 AMOD(PI,VEL) = DMOD(PI,VEL) IF(ILIST .GE. 1) WRITE(IPR,1000) UNB00330 ZERO = 0D0 X60 = 60D0 ISAT = SAT IELEV = ECA IAZ = DCA IQUAD = QUAD IQUAD = KWAD(IQUAD) INJDAY = DAYINJ INJHR = TINJ / 60. INJMIN = AMOD(TINJ,X60) LODAY = DLOCK LOHR = TLOCK/60. LOMIN = AMOD(TLOCK,X60) RSIGO = DOPSIG SIGO = RSIGO**2 NITER = XITER NIT = 2 NIT = 1 NDOP = ELEMD COCODE = 1. MINDOP = DOPMIN C COMPUTE APRIORI WEIGHT MATRIX WT CALL PROP(COVGP,GPROT,SXYZ,3,3) DO 4 I = 1,3 DO 4 J = 1,3 WT(I,J) = 0D0 IF(I.LT.3 .AND. J.LT.3) WT(I,J) = COVGP(I,J) 4 CONTINUE WT(3,3) = FGOSIG**2 IF(ILIST .GE. 3) WRITE(IPR,1023) WT CALL CHOLD(WT,3,3,DETA,ICHK) IF(ICHK .NE. 0) GO TO 190 IF(ILIST .GE. 3) WRITE(IPR,1021) WT IF(ILIST .GE. 2) WRITE(IPR,1010) UNB00390 C INITIALIZATION UNB00340 DO 5 K = 1,3 UNB00370 5 X(K) = 0. UNB00380 DF = 0D0 C ITERATION LOOP UNB00400 DO 80 ITER = 1,NITER UNB00410 C COMPUTE DESIGN MATRIX BASE AND KNOWN FUNCTION F UNB00420 NOBS = 0 UNB00430 T = 0. IF(TCODE .EQ. 2.) T = T + DELAY * 1D-6 WAVE = FO * (1D0 + (FG + DF) * 1D-10) / VEL CALL SATXYZ(T,XS,S1,4) DSDX1 = -(XS(1) - X(1)) / S1 DSDY1 = -(XS(2) - X(2)) / S1 I = 0 IF(ILIST.GE.2.AND.ITER.LE.NIT)WRITE(IPR,1009)I,S1,DSDX1,DSDY1 DO 70 I = 1,NDOP T = T + DOP3(I) CALL SATXYZ(T,XS,S2,4) DSDX2 = -(XS(1) - X(1)) / S2 DSDY2 = -(XS(2) - X(2)) / S2 IF(DOP8(I) .GT. 0.) GO TO 17 IF(ILIST.GE.2.AND.ITER.LE.NIT)WRITE(IPR,1009)I,S2,DSDX2,DSDY2 GO TO 60 C REJECT MISCLOSURES GT 'YMISCL' UNB00600 17 CONTINUE XMISCL = DOP4(I)-DOP3(I)*(OFFREQ+FG+DF)*FO/1D10-WAVE*(S2-S1) IF( ABS(XMISCL) .LT. YMISCL) GO TO 20 UNB00630 IF(ILIST.GE.2.AND.ITER.LE.NIT)WRITE(IPR,1013)I,S2,DSDX2,DSDY2, $ DOP4(I),XMISCL DOP8(I) = -10. GO TO 60 UNB00660 20 NOBS = NOBS + 1 UNB00670 F(NOBS) = XMISCL UNB00680 BASE(1,NOBS) = WAVE * (DSDX2 - DSDX1) BASE(2,NOBS) = WAVE * (DSDY2 - DSDY1) BASE(3,NOBS) = (DOP3(I) + (S2 - S1) / VEL) * FO / 1D10 IF(ILIST.GE.2.AND.ITER.LE.NIT)WRITE(IPR,1008)I,S2,DSDX2,DSDY2, $ NOBS,DOP4(I),F(NOBS),(BASE(K,NOBS),K=1,3) 60 S1 = S2 DSDX1 = DSDX2 DSDY1 = DSDY2 70 CONTINUE UNB00770 C LEAST SQUARES APPROXIMATION UNB00780 IF(NOBS .LT. MINDOP) GO TO 140 UNB00790 CALL LSA(F,NOBS,BASE,3,3,1,1,W,C,R,SIGMA,COV,ICHK) IF(ICHK .NE. 0) GO TO 145 X(1) = X(1) + C(1) X(2) = X(2) + C(2) DF = DF + C(3) DO 71 I = 1,3 71 SD(I) = SQRT(COV(I,I)) STDDEV = SQRT(SIGMA) UNB00860 ALF = 1D0 - (1D0 - ALF1) / NOBS UNB00870 RTEST = XNORM(ALF,ZERO,SQRT(DOPSIG)) IF(ILIST .GE. 2) WRITE(IPR,1015) ALF,RTEST UNB00890 IF(ILIST.GE.1) WRITE(IPR,1003) ITER,STDDEV,NOBS,X(1),X(2),DF,SD IF(ILIST .GE. 2) WRITE(IPR,1007) (R(I),I=1,NOBS) UNB00910 C UNPACK RESIDUALS TO COMPARE WITH DOPPLER ARRAY UNB00920 NRES = 0 NREJ = 0 DO 75 I = 1,NDOP DOP9(I) = 0. IF(DOP8(I) .LE. 0.) GO TO 75 NRES = NRES + 1 IF(NRES .GT. NOBS) GO TO 130 DOP9(I) = R(NRES) IF(ABS(DOP9(I)) .LE. RTEST) GO TO 75 DOP8(I) = -9. NREJ = NREJ + 1 75 CONTINUE C CONVERGENCE TEST UNB01090 IF( ABS(C(1)) .LT. EPS1 .AND. ABS(C(2)) .LT. EPS1 .AND. UNB01100 * ABS(C(3)) .LT. EPS2 .AND. NREJ .EQ. 0) GO TO 90 UNB01110 80 CONTINUE UNB01120 C REJECT IF TOO MANY ITERATIONS UNB01130 GO TO 150 UNB01140 C WRITE DOPPLERS AND RESIDUALS UNB01150 90 CONTINUE UNB01160 IF(ILIST .GE. 1) WRITE(IPR,1024) BUG***** IF(ILIST.GE.1) WRITE(IPR,1004) (I,DOP3(I),DOP4(I),DOP7(I), $ DOP8(I),DOP9(I),I=1,NDOP) C REJECT PASS IF SIGMA FAILS CHI-SQUARE TEST AT 100*ALF2% CONFIDENCE LEVUNB01200 STEST = CHISQ(ALF2,NOBS-3) * SIGO / (NOBS-3) UNB01210 IF(SIGMA .GT. STEST) GO TO 170 UNB01220 C REJECT IF FINAL X,Y TOO LARGE UNB01230 IF( ABS(X(1)) .GT. EPS3 .OR. ABS(X(2)) .GT. EPS3) GO TO 160 UNB01240 C REJECT IF STD DEVS OF SOLUTION TOO LARGE UNB01250 IF(SD(1).GT.EPS4 .OR. SD(2).GT.EPS4 .OR. SD(3).GT.EPS5) GO TO 180 UNB01260 C STORE ACCEPTED PASSES UNB01270 IF(NULL .NE. 0) RETURN UNB01280 IKEPT = IKEPT + 1 UNB01290 C PRINT ONE LINE SUMMARY OF PASS UNB01300 WRITE(IPR,1002) NPASS,IKEPT,ITER,SIGMA,ISAT,IELEV,IAZ,IQUAD, $ INJDAY,INJHR,INJMIN,LODAY,LOHR,LOMIN,NDOP,NOBS,X(1),X(2),DF,SD PHEAD(1) = IKEPT BUG***** I = 56 PHEAD(I) = ITER PHEAD(I+1) = NOBS PHEAD(I+2) = SIGMA PHEAD(I+3) = X(1) PHEAD(I+4) = X(2) PHEAD(I+5) = DF PHEAD(I+6) = SD(1) PHEAD(I+7) = SD(2) PHEAD(I+8) = SD(3) C ROTATE SOLUTION INTO AVERAGE TERRESTRIAL COORDINATES IF(ILIST .GE. 1) WRITE(IPR,1020) COVGP DO 105 I = 1,2 DO 105 J = 1,2 105 COVGP(I,J) = COV(I,J) IF(ILIST .GE. 1) WRITE(IPR,1022) COVGP DO 110 I = 1,3 XAT(I) = 0D0 DO 110 J = 1,3 XAT(I) = XAT(I) + GPROT(J,I) * X(J) COVAT(I,J) = 0D0 DO 110 K = 1,3 DO 110 L = 1,3 110 COVAT(I,J) = COVAT(I,J) + GPROT(K,I)*GPROT(L,J)*COVGP(K,L) IF(ILIST .GE. 1) WRITE(IPR,1017) XAT,COVAT C ROTATE SOLUTION INTO GEODETIC COORDINATES DO 120 I = 1,3 XPLH(I) = 0D0 DO 120 J = 1,3 XPLH(I) = XPLH(I) + TOPROT(I,J) * XAT(J) COVPL(I,J) = 0D0 DO 120 K = 1,3 DO 120 L = 1,3 120 COVPL(I,J) = COVPL(I,J)+TOPROT(I,K)*TOPROT(J,L)*COVAT(K,L) IF(ILIST .GE. 1) WRITE(IPR,1018) XPLH,COVPL RETURN UNB01450 C ERROR RETURNS UNB01460 130 WRITE(IPR,1005) UNB01470 STOP UNB01480 140 WRITE(IPR,1006) NPASS,MINDOP,ISAT,IELEV,IAZ,IQUAD,INJDAY,INJHR, UNB01490 $ INJMIN,LODAY,LOHR,LOMIN,NDOP,NOBS NULL = 2 RETURN UNB01520 145 NULL = 3 RETURN 150 WRITE(IPR,1011) NPASS,NITER,ISAT,IELEV,IAZ,IQUAD,INJDAY,INJHR, UNB01530 $ INJMIN,LODAY,LOHR,LOMIN,NDOP,NOBS NULL = 3 RETURN UNB01560 160 WRITE(IPR,1012) NPASS,ISAT,IELEV,IAZ,IQUAD,INJDAY,INJHR, $ INJMIN,LODAY,LOHR,LOMIN,NDOP,NOBS NULL = 3 RETURN UNB01600 170 WRITE(IPR,1014) NPASS,SIGMA,ISAT,IELEV,IAZ,IQUAD,INJDAY,INJHR, UNB01610 $ INJMIN,LODAY,LOHR,LOMIN,NDOP,NOBS NULL = 3 RETURN UNB01640 180 WRITE(IPR,1016) NPASS,ISAT,IELEV,IAZ,IQUAD,INJDAY,INJHR, $ INJMIN,LODAY,LOHR,LOMIN,NDOP,NOBS NULL = 3 RETURN UNB01680 190 WRITE(IPR,1019) STOP C FORMAT STATEMENTS UNB01690 1000 FORMAT(1X,30HSUBROUTINE GUIER ENTERED ) UNB01700 1002 FORMAT(1X,2I4,I3,F5.1,2I3,I5,1X,A4,1X,2(I4,2I3),2X,2I4,6F7.1) BUG***** 1003 FORMAT(9X,I3,F10.5,10X,I4,3F10.3,10X,3F10.3) 1004 FORMAT(5X,I5,F10.6,F13.3,F6.1,F6.0,F8.3) BUG***** 1005 FORMAT(10X,25HMOBS GT NOBS IN GUIER ) UNB01740 1006 FORMAT(I5,5X, 9HREJECT-LT,I4,7H DOPS ,3I3,2X,A4,2X,2(I4,2I3),2X, UNB01750 * 2I4) UNB01760 1007 FORMAT(5X,10F11.3) 1008 FORMAT(5X,I3,F11.1,2F11.6,I3,F13.3,F8.3,2F11.6,F11.2) 1009 FORMAT(5X,I3,F11.1,2F11.6) 1010 FORMAT(4X,7HMSG ROW,5X,1HS,8X,4HDSDX,7X,4HDSDY,4X,12HNOBS DOPVAC,UNB01800 * 5X,1HF,15X,10HBASE 1 2 3) UNB01810 1011 FORMAT(I5,5X,9HREJECT-GT,I4,7H ITERS ,3I3,2X,A4,2X,2(I4,2I3),2X, UNB01820 * 2I4) UNB01830 1012 FORMAT(I5,5X,16HREJECT-FIX NE IN,4X,3I3,2X,A4,2X,2(I4,2I3),2X,2I4)UNB01840 1013 FORMAT(5X,I3,F11.1,2F11.6,3X,F13.3,E13.4,21H MISCLOSURE TOO LARGE) 1014 FORMAT(I5,5X,12HREJECT-CHISQ ,F8.3,3I3,2X,A4,2X,2(I4,2I3),2X,2I4) 1015 FORMAT(5X,9HALF/RTEST ,2F10.3) 1016 FORMAT(I5,5X,16HREJECT-HI STDDEV,4X,3I3,2X,A4,2X,2(I4,2I3),2X,2I4)UNB01880 1017 FORMAT(5X,30HAVG TERR COORD CORRECTIONS ,3F13.3/3(/35X,3F13.3)) 1018 FORMAT(5X,30HGEODETIC COORD CORRECTIONS (M),3F13.3/3(/35X,3F13.3)) 1019 FORMAT(5X,30HWT MATRIX IN GUIER SINGULAR ) 1020 FORMAT(5X,30HAPRIORI GUIER PLANE COORD COV ,3(/35X,3F13.3)) 1021 FORMAT(5X,30HGUIER PLANE WEIGHT MATRIX ,3(/35X,3F13.3)) 1022 FORMAT(5X,30HEST GUIER PLANE COORD COV ,3(/35X,3F13.3)) 1023 FORMAT(5X,30HGUIER PLANE COV MATRIX ,3(/35X,3F13.3)) 1024 FORMAT(9X,6HI DOP3,11X,22HDOP4 DOP7 DOP8 DOP9 ) BUG***** END UNB01890 SUBROUTINE SPASS C STORE PASS HEADER, ORBIT ARRAYS, AND DOP1 AND DOP2 ARRAYS C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), $ DOP7(260),DOP8(260),DOP9(260) EQUIVALENCE (PHEAD(26),DOPNUM) IF(ILIST .GE. 1) WRITE(IPR,1000) NDOP = DOPNUM IF(NULL .NE. 0) GO TO 5 WRITE(IPN,1001) PHEAD WRITE(IPN,1002) ORB1 WRITE(IPN,1004) ORB3 NDOP1 = NDOP + 1 WRITE(IPN,1008) (I,DOP1(I),DOP2(I),DOP3(I), $ DOP4(I),DOP5(I),DOP6(I),DOP7(I),DOP8(I),DOP9(I),I=1,NDOP1) IF(ILIST .GE. 3) WRITE(IPR,1001) PHEAD IF(ILIST .GE. 3) WRITE(IPR,1002) ORB1 IF(ILIST .GE. 3) WRITE(IPR,1003) ORB2 IF(ILIST .GE. 3) WRITE(IPR,1004) ORB3 IF(ILIST.GE.3)WRITE(IPR,1008) (I,DOP1(I),DOP2(I),DOP3(I), $ DOP4(I),DOP5(I),DOP6(I),DOP7(I),DOP8(I),DOP9(I),I=1,NDOP1) RETURN 5 CONTINUE RETURN C ERROR RETURNS 200 CONTINUE RETURN C FORMAT STATEMENTS 1000 FORMAT(1X,30HSUBROUTINE SPASS ENTERED ) 1001 FORMAT(1X,F5.0,F4.0,F5.0,F6.0,F5.0,F6.0,F12.6,F5.1,2F5.0,F3.0, $ F6.1,F7.1,F5.1/1X,2F5.1,6F4.0,6F5.0,F3.0/1X,3F13.3/1X,9F8.0/ $ 1X,6F13.10/1X,3F13.3,F4.0,F3.0,F12.6,F7.3,F4.0/ $ 1X,F3.0,F5.0,F5.1,6F7.1/1X,36F2.0) 1002 FORMAT(1X,5F15.7) 1003 FORMAT(1X,11F7.0) 1004 FORMAT(1X,8F9.4) 1007 FORMAT(1X,3(F10.6,F11.3,F5.0)) 1008 FORMAT(1X,I3,2F11.2,F10.6,F12.3,F7.3,F8.3,F5.1,F5.0,F7.3) END SUBROUTINE ROTREF(NUM,NAXIS,ANGLE,ROT) UNB00020 C COMPUTE PRODUCT MATRIX OF A SEQUENCE OF ROTATIONS AND REFLECTIONS UNB00030 C INPUT ARGUMENTS UNB00040 C NUM = NUMBER OF ROTATIONS AND REFLECTIONS IN SEQUENCE UNB00050 C NAXIS(NUM) = SEQUENCE OF ROTATION AND REFLECTION AXES UNB00060 C FOR ROTATIONS USE 1,2, OR 3 UNB00070 C FOR REFLECTIONS USE -1, -2, OR -3 UNB00080 C ANGLE(NUM) = SEQUENCE OF ROTATION ANGLES IN RADIANS UNB00090 C FOR REFLECTIONS THIS ANGLE IS IGNORED (SET TO ZERO) UNB00100 C OUTPUT ARGUMENT UNB00110 C ROT(3,3) = PRODUCT MATRIX UNB00120 DIMENSION ROT(3,3),R1(3,3),R2(3),ANGLE(NUM),NAXIS(NUM) UNB00130 C**MACHINE DEPENDENT STATEMENTS UNB00140 DOUBLE PRECISION ROT,R1,R2,ANGLE,EPS,COS,SIN,ABS,DCOS,DSIN,DABS UNB00150 COS(EPS) = DCOS(EPS) UNB00160 SIN(EPS) = DSIN(EPS) UNB00170 ABS(EPS) = DABS(EPS) UNB00180 EPS = 1D-15 UNB00190 C SET ROT = IDENTITY MATRIX UNB00200 DO 1 I = 1,3 UNB00210 DO 1 J = 1,3 UNB00220 ROT(I,J) = 0. UNB00230 IF(I .EQ. J) ROT(I,J) = 1. UNB00240 1 CONTINUE UNB00250 C CHECK ELEMENTS OF NAXIS AND SET REFLECTION ELEMENTS OF ANGLE = 0 UNB00260 DO 2 N = 1,NUM UNB00270 IF(NAXIS(N) .EQ. 0 .OR. NAXIS(N) .LT. -3 .OR. NAXIS(N) .GT. 3) UNB00280 * GO TO 5 UNB00290 IF(NAXIS(N) .LT. 0) ANGLE(N) = 0. UNB00300 2 CONTINUE UNB00310 C PROCESS SEQUENCE OF ROTATIONS AND REFLECTIONS ONE AT A TIME UNB00320 DO 4 N = 1,NUM UNB00330 C DEFINE THREE AXES FOR CURRENT ROTATION OR REFLECTION UNB00340 N1 = IABS(NAXIS(N)) UNB00350 N2 = MOD(N1,3) + 1 UNB00360 N3 = MOD(N2,3) + 1 UNB00370 C DEFINE DIAGONAL ELEMENTS UNB00380 R1(N1,N1) = 1. UNB00390 IF(NAXIS(N) .LT. 0.) R1(N1,N1) = -1. UNB00400 R1(N2,N2) = COS(ANGLE(N)) UNB00410 R1(N3,N3) = R1(N2,N2) UNB00420 C DEFINE NON-ZERO OFF-DIAGONAL ELEMENTS UNB00430 R1(N2,N3) = SIN(ANGLE(N)) UNB00440 R1(N3,N2) = - R1(N2,N3) UNB00450 C DEFINE ZERO OFF-DIAGONAL ELEMENTS UNB00460 R1(N1,N2) = 0. UNB00470 R1(N1,N3) = 0. UNB00480 R1(N2,N1) = 0. UNB00490 R1(N3,N1) = 0. UNB00500 C FORM PRODUCT ( SET ROT = R1 * ROT) UNB00510 DO 4 J = 1,3 UNB00520 DO 3 I = 1,3 UNB00530 R2(I) = 0. UNB00540 DO 3 K = 1,3 UNB00550 3 R2(I) = R2(I) + R1(I,K) * ROT(K,J) UNB00560 DO 4 I = 1,3 UNB00570 ROT(I,J) = R2(I) UNB00580 IF(ABS(ROT(I,J)) .LT. EPS) ROT(I,J) = 0. UNB00590 4 CONTINUE UNB00600 RETURN UNB00610 C ERROR RETURN UNB00620 5 WRITE(6,6) N,NAXIS(N) UNB00630 RETURN UNB00640 C FORMAT STATEMENTS UNB00650 6 FORMAT(21H ILLEGAL VALUE NAXIS( ,I3,4H) = ,I5,10H IN ROTREF ) UNB00660 END UNB00670 C UNB00010 SUBROUTINE ROTDOT(NUM,NAXIS,ANGLE,DANGLE,DROT) UNB00020 C COMPUTE TIME DERIVATIVE OF PRODUCT ROTATION MATRIX UNB00030 C INPUT ARGUMENTS UNB00040 C NUM = NUMBER OF ROTATIONS AND REFLECTION IN PRODUCT SEQUENCE UNB00050 C NAXIS(NUM) = SEQUENCE OF ROTATION AND REFLECTION AXES UNB00060 C ROTATION AXES USE 1,2,3 UNB00070 C REFLECTION AXES USE -1,-2,-3 UNB00080 C ANGLE(NUM) = SEQUENCE OF ROTATION ANGLES IN RADIANS UNB00090 C FOR REFLECTIONS SUBROUTINE SETS ANGLE TO ZERO UNB00100 C DANGLE(NUM) = SEQUENCE OF TIME DERIVATIVES OF ELEMENTS OF ANGLE UNB00110 C OUTPUT ARGUMENT UNB00120 C DROT(3,3) = TIME DERIVATIVE OF PRODUCT MATRIX UNB00130 C THIS SUBROUTINE CALLS SUBROUTINE ROTREF UNB00140 C**MACHINE DEPENDENT STATEMENT UNB00150 DOUBLE PRECISION DROT,R1,R2,R3,ANGLE,DANGLE,PI2,ANGLEN UNB00160 DIMENSION DROT(3,3),R1(3,3),R2(3,3),R3(3,3),ANGLE(NUM),DANGLE(NUM)UNB00170 * ,NAXIS(NUM) UNB00180 PI2 = 3.141592653589793D0 / 2D0 UNB00190 C INITIALIZE (SET DROT = ZERO MATRIX) UNB00200 DO 10 I = 1,3 UNB00210 DO 10 J = 1,3 UNB00220 10 DROT(I,J) = 0. UNB00230 C COMPUTE DROT AS SUM OF NUM MATRICES, EACH OF WHICH IS A PRODUCT OF NUMUNB00240 C MATRICES UNB00250 DO 30 N = 1,NUM UNB00260 IF(DANGLE(N) .EQ. 0.) GO TO 30 UNB00270 IF(NAXIS(N) .LE. 0) GO TO 30 UNB00280 ANGLEN = ANGLE(N) UNB00290 NAXISN = NAXIS(N) UNB00300 ANGLE(N) = ANGLEN + PI2 UNB00310 CALL ROTREF(NUM,NAXIS,ANGLE,R1) UNB00320 NAXIS(N) = - (MOD(NAXISN,3) + 1) UNB00330 CALL ROTREF(NUM,NAXIS,ANGLE,R2) UNB00340 NAXIS(N) = -(MOD(NAXISN+1,3) + 1) UNB00350 CALL ROTREF(NUM,NAXIS,ANGLE,R3) UNB00360 DO 20 I = 1,3 UNB00370 DO 20 J = 1,3 UNB00380 20 DROT(I,J) = DROT(I,J) + DANGLE(N)*(R1(I,J) - 0.5*(R2(I,J) + UNB00390 * R3(I,J))) UNB00400 ANGLE(N) = ANGLEN UNB00410 NAXIS(N) = NAXISN UNB00420 30 CONTINUE UNB00430 RETURN UNB00440 END UNB00450 C UNB00010 FUNCTION CHISQ(A,N) UNB00020 C PERCENTILES OF THE CHI-SQUARE DISTRIBUTION X(N) UNB00030 C INPUT UNB00040 C A = PROBABILITY INTEGRAL FROM ZERO TO CHISQ UNB00050 C N = NUMBER OF DEGREES OF FREEDOM UNB00060 C OUTPUT UNB00070 C CHISQ = ABSCISSA VALUE OF X(N) CORRESPONDING TO PROBABILITY A UNB00080 C ABRAMOWITZ AND STEGUN EQN 26.4.17 IS USED (WILSON-HILFERTY APPROX) UNB00090 C ACCURACY BETTER THAN 0.04 FOR N .GT. 1 UNB00100 C**MACHINE DEPENDENT STATEMENTS UNB00110 DOUBLE PRECISION A,CHISQ,XNORM,SQRT,DSQRT UNB00120 SQRT(A) = DSQRT(A) UNB00130 CHISQ =N*(1-2./(9.*N)+XNORM(A,0D0,1D0)* SQRT(2D0/(9.*N)))**3 UNB00140 RETURN UNB00150 END UNB00160 C UNB00010 FUNCTION XNORM(ALF,XMEAN,SIG) UNB00020 C PERCENTILES OF THE NORMAL DISTRIBUTION N(XMEAN,SIG) UNB00030 C INPUT UNB00040 C ALF = PROBABILITY INTEGRAL FROM NEGATIVE INFINITY TO XNORM UNB00050 C XMEAN = POPULATION MEAN UNB00060 C SIG = POPULATION STANDARD DEVIATION UNB00070 C OUTPUT UNB00080 C XNORM = ABSCISSA VALUE OF N(XMEAN,SIG) CORRESPONDING TO UNB00090 C PROBABILITY ALF UNB00100 C ABRAMOWITZ AND STEGUN EQUATION 26.2.23 IS USED UNB00110 C ACCURACY BETTER THAN THAN 0.00045 UNB00120 DIMENSION C(6) UNB00130 C**MACHINE DEPENDENT STATEMENTS UNB00140 DOUBLE PRECISION XNORM,ALF,XMEAN,SIG,P,T,T2,T3,XP,SQRT,ALOG,DSQRT,UNB00150 * DLOG UNB00160 SQRT(P) = DSQRT(P) UNB00170 ALOG(P) = DLOG(P) UNB00180 DATA C(1),C(2),C(3),C(4),C(5),C(6)/2.515517,.802853,.010328, UNB00190 * 1.432788,.189269,.001308/ UNB00200 IF(ALF .GE. 1. .OR. ALF .LE. 0.) GO TO 20 UNB00210 SIGN = -1. UNB00220 P = ALF UNB00230 IF(ALF .LT. .5) GO TO 10 UNB00240 SIGN = 1. UNB00250 P = 1. - ALF UNB00260 10 T = SQRT(ALOG(1./P**2)) UNB00270 T2 = T * T UNB00280 T3 = T * T2 UNB00290 XP = T - (C(1)+C(2)*T+C(3)*T2) / (1.+C(4)*T+C(5)*T2+C(6)*T3) UNB00300 XP = XP * SIGN UNB00310 XNORM = XMEAN + XP * SIG UNB00320 RETURN UNB00330 C ERROR RETURN UNB00340 20 WRITE(6,1001) ALF UNB00350 XNORM = 0. UNB00360 RETURN UNB00370 C FORMAT STATEMENTS UNB00380 1001 FORMAT(25H XNORM INPUT PROBABILITY= ,E20.10) UNB00390 END UNB00400 C UNB00010 C UNB00010 SUBROUTINE LSA(F,N,BASE,M,IRBASE,IDENT,IRW,W,C,R,RNORM,A,ICHK) C LEAST SQUARES APPROXIMATION UNB00020 C INPUT ARGUMENTS UNB00030 C F(N) = INPUT FUNCTION TO BE APPROXIMATED UNB00040 C BASE(M,N) = (IF DIMENSIONED) COMPUTED DESIGN MATRIX UNB00050 C (IF NOT DIMENSIONED) FUNCTION SUBROUTINE NAME TO UNB00060 C COMPUTE BASE FUNCTIONS UNB00070 C IRBASE = ROW DIMENSION OF ARRAY BASE IN CALLING ROUTINE UNB00080 C IDENT = WEIGHT MATRIX SWITCH(IDENT=1 IF WEIGHT MATRIX IS IDENTITY, UNB00090 C IDENT = N OTHERWISE) UNB00100 C IRW = ROW DIMENSION OF ARRAY W IN CALLING ROUTINE ( .EQ. 1 IF UNDIM) UNB00110 C W(N,N) = WEIGHT MATRIX (IGNORED IF IDENT = 1) UNB00120 C OUTPUT ARGUMENTS UNB00130 C C(M) = VECTOR OF LEAST SQUARES APPROXIMATION COEFFICIENTS UNB00140 C R(N) = RESIDUAL VECTOR = F(N) - P(N) UNB00150 C RNORM = VARIANCE FACTOR C A(M,M) = COVARIANCE MATRIX OF C C ICHK = 0 SUCCESSFUL RETURN UNB00180 C 1 SINGULAR (N .LT. M) UNB00190 C 2 A AND B DIMENSIONED .LT. M UNB00200 IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(N),C(M),R(N),A(M,M),W(IRW,IDENT) C**OPTIONAL STATEMENT (SEE COMMENTS ABOVE) UNB00220 DIMENSION BASE(IRBASE,N) UNB00230 C**MACHINE DEPENDENT STATEMENT UNB00240 DOUBLE PRECISION F,BASE,W,C,R,RNORM,SD,A,B,AA,DETA,P,SQRT,DSQRT UNB00250 C**CHOOSE IRDA .GE. M UNB00260 DIMENSION B(100) DATA IRDA/100/ C**MACHINE DEPENDENT STATEMENT UNB00290 SQRT(AA) = DSQRT(AA) UNB00300 ICHK = 0 UNB00310 IF(IRDA .LT. M) GO TO 130 UNB00320 IF(N .LT. M) GO TO 140 UNB00330 C CLEAR ARRAYS A AND B UNB00340 DO 10 I = 1,M UNB00350 B(I) = 0. UNB00360 DO 10 J = 1,M UNB00370 10 A(I,J) = 0. UNB00380 C CONSTRUCT NORMAL EQUATION ARRAYS A AND B UNB00390 DO 60 IT1 = 1,N UNB00400 DO 60 J = 1,M UNB00410 AA = BASE(J,IT1) UNB00420 IF(AA .EQ. 0.) GO TO 60 UNB00430 IF(IDENT .NE. 1) GO TO 30 UNB00440 B(J) = B(J) + AA * F(IT1) UNB00450 DO 20 K = J,M UNB00460 20 A(K,J) = A(K,J) + AA * BASE(K,IT1) UNB00470 GO TO 60 UNB00480 30 DO 50 IT2 = 1,N UNB00490 IF(W(IT1,IT2) .EQ. 0.) GO TO 50 UNB00500 B(J) = B(J) + AA * F(IT2) * W(IT1,IT2) UNB00510 DO 40 K = J,M UNB00520 40 A(K,J) = A(K,J) + AA * BASE(K,IT2) * W(IT1,IT2) UNB00530 50 CONTINUE UNB00540 60 CONTINUE UNB00550 C CHOLESKI INVERSION OF A UNB00560 CALL CHOLD(A,M,M,DETA,JCHK) IF(JCHK .NE. 0) GO TO 200 UNB00580 C COMPUTE COEFFICIENTS C UNB00590 DO 70 I = 1,M UNB00600 C(I) = 0. UNB00610 DO 70 J = 1,M UNB00620 70 C(I) = C(I) + A(I,J) * B(J) UNB00630 C COMPUTE APPROXIMANT P AND RESIDUAL R UNB00640 DO 90 I = 1,N UNB00650 P = 0. UNB00660 DO 80 J = 1,M UNB00670 80 P = P + C(J) * BASE(J,I) UNB00680 90 R(I) = F(I) - P UNB00690 C COMPUTE QUADRATIC NORM RNORM OF R UNB00700 RNORM = 0. UNB00710 IF(IDENT .NE. 1) GO TO 110 UNB00720 DO 100 IT1 = 1,N UNB00730 100 RNORM = RNORM + R(IT1) * R(IT1) UNB00740 GO TO 121 UNB00750 110 DO 120 IT1 = 1,N UNB00760 DO 120 IT2 = 1,N UNB00770 IF(W(IT1,IT2) .NE. 0.) RNORM = RNORM + R(IT1) * R(IT2) * UNB00780 * W(IT1,IT2) UNB00790 120 CONTINUE UNB00800 121 CONTINUE UNB00810 RNORM = RNORM / (N-M) C COMPUTE COVARIANCE MATRIX OF COEFFICIENTS DO 125 I = 1,M DO 125 J = 1,M 125 A(I,J) = A(I,J) * RNORM RETURN UNB00850 C ERROR RETURN UNB00860 130 WRITE(6,1001) UNB00870 ICHK = 2 UNB00880 RETURN UNB00890 140 WRITE(6,1002) UNB00900 200 ICHK = 1 UNB00910 RETURN UNB00920 C FORMAT STATEMENTS UNB00930 1001 FORMAT(50H ARRAYS A AND B IN LSA DIMENSIONED .LT. ARGUMENT M ) UNB00940 1002 FORMAT(54H NUMBER OF VALUES F(N) .LT. NUMBER OF COEF C(M) IN LSA ,UNB00950 * 18H - SINGULAR SYSTEM ) UNB00960 END UNB00970 SUBROUTINE CHOLD(A,IRDA,NA,DETA,JCHK) UNB00020 C MATRIX INVERSION USING CHOLESKI DECOMPOSITION UNB00030 C INPUT ARGUMENTS UNB00040 C A = ARRAY CONTAINING POSITIVE DEFINITE SYMMETRIC INPUT MATRIX UNB00050 C IRDA = ROW DIMENSION OF ARRAY CONTAINING INPUT MATRIX UNB00060 C NA = SIZE OF INPUT MATRIX UNB00070 C OUTPUT ARGUMENTS UNB00080 C A = CONTAINS INVERSE OF INPUT MATRIX (INPUT DESTROYED) UNB00090 C DETA = DETERMINANT OF INPUT MATRIX UNB00100 C JCHK = 0 SUCCESSFUL RETURN UNB00110 C 1 SINGULAR MATRIX UNB00120 C 2 MATRIX OF ZERO DIMENSION UNB00130 DIMENSION A(IRDA,NA) UNB00140 C**MACHINE DEPENDENT STATEMENTS UNB00150 DOUBLE PRECISION A,DETA,SUM,SQRT,ABS,SING,DSQRT,DABS UNB00160 SQRT(DETA) = DSQRT(DETA) UNB00170 ABS(DETA) = DABS(DETA) UNB00180 SING = 1D-20 UNB00190 JCHK = 0 UNB00200 C CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX UNB00210 IF(NA .LT. 1) GO TO 18 UNB00220 DETA = A(1,1) UNB00230 A(1,1) = SQRT(A(1,1)) UNB00240 IF(NA .EQ. 1) GO TO 6 UNB00250 DO 1 I = 2,NA UNB00260 1 A(I,1) = A(I,1) / A(1,1) UNB00270 DO 5 J = 2,NA UNB00280 SUM = 0. UNB00290 J1 = J - 1 UNB00300 DO 2 K = 1,J1 UNB00310 2 SUM = SUM + A(J,K) ** 2 UNB00320 DETA = DETA * (A(J,J) - SUM) UNB00330 A(J,J) = SQRT(A(J,J) - SUM) UNB00340 IF(J .EQ. NA) GO TO 5 UNB00350 J2 = J + 1 UNB00360 DO 4 I = J2,NA UNB00370 SUM = 0. UNB00380 DO 3 K = 1,J1 UNB00390 3 SUM = SUM + A(I,K) * A(J,K) UNB00400 4 A(I,J) = (A(I,J) - SUM) / A(J,J) UNB00410 5 CONTINUE UNB00420 6 IF(ABS(DETA) .LT. SING) GO TO 16 UNB00430 C INVERSION OF LOWER TRIANGULAR MATRIX UNB00440 DO 7 I = 1,NA UNB00450 7 A(I,I) = 1. / A(I,I) UNB00460 IF(NA .EQ. 1) GO TO 10 UNB00470 N1 = NA - 1 UNB00480 DO 9 J = 1,N1 UNB00490 J2 = J + 1 UNB00500 DO 9 I = J2,NA UNB00510 SUM = 0. UNB00520 I1 = I - 1 UNB00530 DO 8 K = J,I1 UNB00540 8 SUM = SUM + A(I,K) * A(K,J) UNB00550 9 A(I,J) = - A(I,I) * SUM UNB00560 C CONSTRUCTION OF INVERSE OF INPUT MATRIX UNB00570 10 DO 15 J = 1,NA UNB00580 IF(J .EQ. 1) GO TO 12 UNB00590 J1 = J - 1 UNB00600 DO 11 I = 1,J1 UNB00610 11 A(I,J) = A(J,I) UNB00620 12 DO 14 I = J,NA UNB00630 SUM = 0. UNB00640 DO 13 K = I,NA UNB00650 13 SUM = SUM + A(K,I) * A(K,J) UNB00660 14 A(I,J) = SUM UNB00670 15 CONTINUE UNB00680 RETURN UNB00690 C ERROR RETURNS UNB00700 16 WRITE(6,17) DETA UNB00710 JCHK = 1 UNB00720 RETURN UNB00730 18 WRITE(6,19) UNB00740 JCHK = 2 UNB00750 RETURN UNB00760 C FORMAT STATEMENTS UNB00770 17 FORMAT(33H SINGULAR MATRIX IN CHOLD. DET = ,E20.5) UNB00780 19 FORMAT(35H MATRIX OF DIMENSION ZERO IN CHOLD ) UNB00790 END UNB00800