C C C C******************************************************************************* C* * C* * C* ***************************** * C* * * * C* * 'P R E-B A T C H'-PROGRAM * * C* * * * C* ***************************** * C* * C* SATELLITE (XYZ) CARTESIAN COORDINATES IN THE AVERAGE TERRESTRIAL SYSTEM * C* ----------------------------------------------------------------------- * C* PASS IDENTIFICATION AND CHARACTERISTICS W.R.T. THE TRACKING STATION * C* ------------------------------------------------------------------- * C* ACCEPTED DOPPLER COUNTS CORRECTED FOR IONO. & TROPO. REFRACTION * C* --------------------------------------------------------------- * C* FROM A SINGLE MOTIONLESS TRACKING STATION * C* ----------------------------------------- * C* * C* PREPARED BY : * C* 'MOHAMED M. NASSAR' / *G.S.* - (DEC. 1971) * C* * C* THIS PROGRAM IS A MODIFIED VERSION FROM THE ORIGINAL PROGRAM AVAILABLE IN * C* THE SURVEYING ENGINEERING DEPARTMENT (UNB) , WHICH IS : * C* * C* UNB DOPPLER SATELLITE PROGRAM ONE - 'PREPROCESSOR' * C* WRITTEN BY D. WELLS (JANUARY 1971) * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* * C* THE PURPOSE OF THIS PROGRAM IS TO PREPARE THE INPUT DATA (ON PUNCHED CARDS) * C* FOR THE OTHER PROGRAM 'BATCH' ,WHICH PERFORMS THE LEAST SQUARES ADJUSTMENT * C* AND SOLVE FOR THE COORDINATES OF THE TRACKING STATION. SPECIFICALLY THE PRE* C* -BATCH PROGRAM PROVIDES THE (X,Y,Z) SATELLITE COORDINATES IN THE AVERAGE * C* TERRESTRIAL SYSTEM ,AND THE CORRECTED DOPPLERS ; THIS IS OBTAINED AFTER EXT-* C* ENSIVE TREATMENT OF THE MAJORITY VOTTED DATA OF THE SATELLITE MESSAGE BEING * C* INPUTTED TO THIS PROGRAM. * C* * C* THIS PROGRAM HAS A SUBROUTINE CALLED 'FREAD' ,WHICH PERMITS A FREE-READING * C* OF ANY INPUTTED NUMERICAL NUMBERS ,WHICH MAKES THE USER HAPPY ABOUT IT ; * C* SINCE THERE IS NO BOTHERING IN FOLLOWING A SPECIAL FORMAT IN PREPARING THE * C* INPUT DATA. UNFORTUNATELY ,THE 'FREAD' ROUTINE IS UNABLE TO READ ANY CHARAC-* C* TERS ,SO WHEN DEALING WITH ANY CHARACTER MANIPULATION ONE HAS TO USE AN * C* ORDINARY FORMATTED-STAEMENT WHENEVER IT IS REQUIRED. IN THIS PROGRAM SUCH * C* FORMATTED STATEMENT IS USED IN TWO-PLACES TO READ THE NAMES OF THE TRACKING * C* STATION ,THE GEODETIC ELLIPSOID ,AND THE GEOCENTRIC ELLIPSOID. RESPECTIVELY * C* * C* THE CHARACTERISTICS OF THE INPUT AND OUTPUT DATA ARE CAREFULLY AND FULLY * C* EXPLAINED WITHIN THE PROGRAM LISTINGS ,IT MIGHT BE WORTHWHILE HERE JUST TO * C* SUMMARIZE THEIR SEQUENCE AS FOLLOWS : * C* * C* THE REQUIRED INPUT DATA IN THE FOLLOWING SEQUENCE : * C* ------------------------------------------------- * C* * C* 1-THE PROGRAM OPTIONS (DTPERG,HORREJ,NCOPYS,IPUNCH,IDCORR) ... 1 CARD * C C C* 2-NAMES OF THE CHOZEN GEODETIC AND THE GEOCENTRIC ELLIPSOIDS ..... 1 CARD * C* ACCORDING TO THE FOLLOWING FORMAT ; * C* COLS 1 TO 3 : GEODETIC ELLIPSOID NAME .............. A3 * C* COLS 6 TO 8 : GEOCENTRIC ELLIPSOID NAME ............. A3 * C* 3-PARAMETERS AND TRANSLATIONS OF THE GEODETIC ELLIPSOID ; (A,FINV,X0,Y0, * C* Z0 & IDTRAN) ................................ 1 CARD * C* 4-PARAMETERS OF THE GEOCENTRIC ELLIPSOID ; (ACG & FINVCG) ,IF ANY....1 CARD * C* 5-TRACKING STATION-NUMBER (ICHOOZ) .......................... 1 CARD * C* 6-TRACKING STATION NAME (STNAME) ,AND THE DIRECTION OF RECKONING THE LONGIT-* C* UDE (LONDIR) ......... 1 CARD ,ACCORDING TO THE FOLLOWING FORMAT ; * C* COLS 1 TO 4 : TRACKING STATION NAME ................. A4 * C* COLS 7 TO 10 : LONGITUDE DIRECTION ................... A4 * C* 7-APPROXIMATE GEOGRAPHICAL COORDINATES OF THE TRACKING STATION ; LATITUDE & * C* LONGITUDE IN (DEG,MIN,SEC) ,AND ELLIPSOIDAL HEIGHT IN METERS ......1 CARD * C* 8-THE YEAR OF OBSERVATIONS (IYEAR) ,AND THE APRIORI VARIANCE-FACTOR (IAVF) * C* .................. 1 CARD * C* * C* FOR EACH PASS THE FOLLOWING INPUT IS REQUIRED ; * C* 1-THE SATELLITE IDENTIFICATION CARD (DAY # ,PASS # ,SAT # & LOCKON HOUR AND * C* MIN) ............................ 1 CARD * C* 2-THE UNSCALED 11 FIXED ORBITAL PARAMETERS ................... 1 CARD * C* 3-THE SETS OF UNSCALED VARIABLE ORBITAL PARAMETERS (UP TO 17 SETS ,EACH IS * C* THREE PARAMETERS CAN BE READ-IN) ......................... 2 CARDS * C* 4-A DELIMITER 999 ................................. 1 CARD * C* 5-THE TRACKING STATION NUMBER (ISTN) ....................... 1 CARD * C* 6-THE PAIRS OF DOPPLER COUNTS ; THE '400-MHZ' AND THE '150-MHZ' RESPECTIVELY* C* ,EITHER ITT OR MAGNAVOX FORMAT IS ACCEPTED. NOTE THAT THE NUMBER OF THESE * C* DOPPLER PAIRS SHOULD BE LESS BY 2 THAN THE NUMBER OF SETS OF VARIABLE * C* ORBITAL PARAMETERS ................................ 2 CARDS * C* 7- A DELIMITER 99 ...................................... 1 CARD * C* * C* FINALLY ,AT THE END OF DATA FOR ALL PASSES ; PUNCH A DELIMITER 9999 * C* ..................... 1 CARD * C* * C* THE PRINTED OUTPUT WILL INCLUDE THE FOLLOWING : * C* --------------------------------------------- * C* * C* I-THE DESIRED OPTIONS ,THE PARAMETERS OF THE REFERENCE ELLIPSOIDS ,AND THE * C* APPROXIMATE COORDINATES OF THE TRACKING STATION. * C* II-FOR EACH PASS BEING PROCESSED ,THE PROGRAM PRINTS THE FOLLOWING : * C* 1-THE UNSCALED AND SCALED (FIXED & VARIABLE) ORBITAL PARAMETERS. * C* 2-SATELLITE CARTESIAN COORDINATES IN ITS ORBIT-PLANE SYSTEM. * C* 3-SATELLITE COORDINATES IN THE AVERAGE TERRESTRIAL SYSTEM. * C* 4-SATELLITE SUBPOINT GEOGRAPHICAL COORDINATES W.R.T. THE GEODETIC ELLIPS. * C* * C*III-REGARDING THE OBSERVED DATA FROM THE TRACKING STATION FOR EACH PASS ; * C* THE PROGRAM PRINTS THE FOLLOWING : * C* 1-THE SATELLITE CARTESIAN COORDINATES IN THE LOCAL GEODETIC SYSTEM. * C* 2-SATELLITE TOPOCENTRIC COORDINATES ; RANGE , AZIMUTH & ELEVATION ANGLE. * C* 3-SUMMARY OF INFORMATION AT THE SATELLITE CLOSEST APPROACH. * C* 4-DETAILS OF ATMOSPHERIC REFRACTION CORRECTIONS TO THE OBSERVED DOPPLERS. * C* 5-WHETHER THE PASS IS ACCEPTED OR REJECTED. * C* 6-THE NUMBER OF AVAILABLE DOPPLER EQUATIONS FOR LEAST SQUARES ADJUSTMENT. * C* * C* IV-AFTER PROCESSING ALL INPUTTED PASSES ,THE PROGRAM PRINTS A SUMMARY OF * C* THE INPUT/OUTPUT DATA (AS MANY COPIES AS REQUIRED) ; CONTAINING :- * C C C* 1-PARAMETERS OF THE CHOZEN REFERENCE ELLIPSOIDS. * C* 2-APPROXIMATE COORDINATES (BOTH GEOGRAPHICAL & CARTESIAN) OF TRACKING STN. * C* 3-IDENTIFICATION AND CHARACTERISTICS OF ALL ACCEPTED PASSES. * C* 4-SATELLITE CARTESIAN COORDINATES IN THE AVERAGE TERRESTRIAL SYSTEM ,WITH * C* THE CORRESPONDING CORRECTED & ACCEPTED DOPPLERS. * C* * C* V-WHEN DESIRED ,THE PROGRAM WILL PRINT AND PUNCH (FOR 'BATCH'-PROGRAM) THE * C* FOLLOWING DATA :- * C* 1-# OF ACCEPTED PASSES ,# OF OBSERVED SATELLITE POSITIONS ,ROW & COLUMN * C* DIMENSIONS OF THE DESIGN MATRIX-A ,YEAR OF OBSERVATIONS ,AND THE APRIORI * C* VARIANCE-FACTOR. * C* 2-THE WEIGHT COEFFICIENT MATRIX OF THE OBSERVED DOPPLERS (EITHER IDENTITY * C* OR WITH CORRELATION). * C* 3-THE VECTOR OF INITIAL APPROXIMATE VALUES OF THE UNKNOWN PARAMETERS. * C* 4-THE SATELLITE CARTESIAN COORDINATES IN THE AVERAGE TERRESTRIAL SYSTEM , * C* AND THE CORRESPONDING DOPPLERS ; FOR EACH ACCEPTED PASS. * C* 5-THE IDENTIFICATION AND CHARACTEISTICS OF EACH ACCEPTED PASS. * C* * C* 6-NAMES AND PARAMETERS OF THE CHOZEN GEODETIC & GEOCENTRIC ELLIPSOIDS. * C* 7-THE TRACKING STATION NAME & NUMBER. * C* * C* NOTE : THIS PROGRAM IS CAPABLE TO HANDLE UP TO 100 PASSES. * C* * C******************************************************************************* C IMPLICIT REAL*8(A-H,O-Z) NASR DOUBLE PRECISION DSQRT,DABS,DSIN,DCOS,DATAN NASR DIMENSION DEK(17),DAK(17),ETAK(17),MINK(17) NASR DIMENSION X(4),ROT(3,3),ANGLE(3),XS(15,4),NAXIS(3),SUB(15,3) NASR DIMENSION TOPROT(3,3),TOPANG(2),NAXTOP(2),R(15,4) NASR DIMENSION RT(15,4),AZ(15),VA(15) NASR DIMENSION DOPC(14),REFR(14),IDOP(14),DOPR(14),VDOP(14) NASR DIMENSION XOBS(3),XSAT(900,3),IEQN(900),DOP(900) NASR DIMENSION IDL(100),ISPE(100) NASR DIMENSION VAMSUM(100),IAZSUM(100),IDRSUM(100),IZQUAD(4) NASR DIMENSION IQUADR(100),NUMXYZ(100) NASR DIMENSION IEQNON(600),COFACT(600) NASR INTEGER CGELIP,ELNAME/'NIL'/,NEGTIV/'WEST'/ NASR C NASR C***********************************************************************NASR C NASR C INSERT THE VALUES OF THE CONSTANTS NEEDED THROUGH THE EXECUTION NASR C OF THE PROGRAM. NASR C NASR C PHYSICAL CONSTANTS NASR C NASR PI=3.141592653589793D0 NASR WE=4.3752695D-03 NASR DEGRA=PI/180.D0 NASR C NASR C RADIO PROPAGATION CONSTANTS NASR C NASR FREQN=400.D6 NASR DELF=80.D-6 NASR FREQT=FREQN*(1.-DELF) NASR VEL=299792500.D0 NASR C C C WAVENO=FREQN/VEL NASR C NASR C TROPOSPHERIC REFRACTION CONSTANTS NASR C NASR XKD=2.31D0 NASR XKW=0.20D0 NASR THD=(2.5*DEGRA)**2 NASR THW=(1.5*DEGRA)**2 NASR C NASR C***********************************************************************NASR C NASR C READ-IN THE PROGRAM OPTIONS :- NASR C NASR C DTPERG : IS THE CORRECTION FOR THE SCALED TIME OF PERIGEE -IN UT MIN NASR C HORREJ : IS THE HORIZON CUT-OFF ELEVATION LIMIT - IN DEGREES OF ARC NASR C NCOPYS : IS THE REQUIRED NUMBER OF COPIES TO BE PRINTED-OUT FOR THE NASR C 'SUMMARY OF THE (PRE-BATCH)-PROGRAM OUTPUT' NASR C IPUNCH : IS AN INDICATOR FOR PUNCHING AND PRINTING THE NECESSARY NASR C DATA FOR THE (MATLAN)-PROGRAM 'BATCH' ; IF IT EQUALS 1 ,THE NASR C PUNCHING & PRINTING WILL OCCUR. NASR C IDCORR : IS AN INDICATOR TO SPECIFY THE DESIRED PUNCHED WEIGHT COEFF-NASR C ICIENT MATRIX OF THE ACCEPTED DOPPLER COUNTS NEEDED FOR THE NASR C ADJUSTMENT ,WHETHER IT IS AN IDENTITY MATRIX -OR A MATRIX NASR C WITH OFF-DIAGONALS DUE TO THE CORRELATION BETWEEN THE NASR C ADJACENT DOPPLERS ; IF IDCORR=1 **CORRELATION EXISTS*** NASR C NASR DTPERG=FREAD(1,ITERM) NASR HORREJ=FREAD(1,ITERM) NASR NCOPYS=FREAD(1,ITERM) NASR IPUNCH=FREAD(1,ITERM) NASR IDCORR=FREAD(1,ITERM) NASR C NASR C PRINT THE DESIRED INPUTTED PROGRAM-OPTIONS. NASR C NASR PRINT2080 NASR 2080 FORMAT('1'//10X,'DESIRED PROGRAM-OPTIONS ARE :'/10X,'-------------NASR 1--------------'/) NASR PRINT2081,DTPERG,HORREJ NASR 2081 FORMAT(5X,'DTPERG =',D9.2,' : IS THE CORRECTION FOR THE SCALED TIMNASR 1E OF PERIGEE -(IN UT MINUTES)'//5X,'HORREJ =',D9.2,' : IS THE HORINASR 2ZON CUT-OFF ELEVATION LIMIT -(IN DEGREES OF ARC)'/) NASR PRINT2082,NCOPYS,IPUNCH NASR 2082 FORMAT(5X,'NCOPYS =',I4,5X,' : IS THE REQUIRED NUMBER OF COPIES FONASR 1R THE PRINTED OUTPUT SUMMARY'//5X,'IPUNCH =',I4,5X,' : IF EQUALS 1NASR 2 ,PUNCH & PRINT THE DATA FOR (BATCH)-PROGRAM WILL OCCUR '/) NASR PRINT2083,IDCORR NASR 2083 FORMAT(5X,'IDCORR =',I4,5X,' : IF EQUALS 1 ,CORRELATION BETWEEN ANASR 1DJACENT DOPPLERS EXISTS '/) NASR C NASR C***********************************************************************NASR C NASR C READ-IN THE NAMES OF THE GEODETIC AND THE GEOCENTRIC ELLIPSOIDS NASR C RESPECTIVELY ; EACH COMPOSED OF THREE ALPHAMERIC CHARACTERS ,START- NASR C ING BY AN ALPHABETIC. NASR C NASR C C C NOTE THAT IF THE GEODETIC ELLIPSOID ONLY WILL BE USED ,THE NAME OF NASR C THE GEOCENTRIC ELLIPSOID SHOULD BE PUNCHED AS ( NIL ) ,TO CONTROL THENASR C EXECUTION. IN THIS CASE WE WILL HAVE ONLY ONE SET OF ELLIPSOIDAL NASR C PARAMETERS TO BE INPUTTED ,WHICH PERTAINS TO THE GEODETIC ELLIPSOID. NASR C NASR READ5001,GELLIP,CGELIP NASR 5001 FORMAT(A3,2X,A3) NASR C NASR C READ-IN THE PARAMETERS OF THE 'GEODETIC' REFERENCE ELLIPSOID ; THE NASR C SEMI-MAJOR AXIS (IN METERS) ,AND THE RECIPROCAL OF FLATTENNING. NASR C NASR A=FREAD(1,ITERM) NASR FINV=FREAD(1,ITERM) NASR C NASR C COMPUTE THE OTHER PARAMETERS OF THE 'GEODETIC' ELLIPSOID ,NAMELY : NASR C ( F , B & E ) NASR C NASR F=1./FINV NASR B=A*(1.-F) NASR E=DSQRT(2.*F-F*F) NASR C NASR C READ-IN THE TRANSLATION-COMPONENTS OF THE 'GEODETIC' REFERENCE ELLIP-NASR C SOID FROM THE GEOCENTRE ,AND THE INDICATOR OF THEIR EXISTANCE IDTRAN NASR C WHICH SHOULD BE PUNCHED AS 1 ;IF ANYONE OF THE THREE-TRANSLATION NASR C COMPONENTS IS NON-ZERO. OTHERWISE IT SHOULD BE ANY OTHER DIGIT. NASR C NASR X0=FREAD(1,ITERM) NASR Y0=FREAD(1,ITERM) NASR Z0=FREAD(1,ITERM) NASR IDTRAN=FREAD(1,ITERM) NASR IF(CGELIP.EQ.ELNAME) GO TO 5094 NASR C NASR C READ-IN THE SEMI -MAJOR AXIS (IN METERS) ,AND RECIPROCAL OF FLATTENN-NASR C ING (DIMENSIONLESS) ; OF THE ARBITRARILY CHOZEN GEOCENTRIC ELLIPSOID NASR C (NOTE THAT THE 3-TRANSLATION COMPONENTS ARE EQUAL TO ZERO ,IN CASE NASR C OF THE 'GEOCENTRIC' ELLIPSOID). NASR C NASR ACG=FREAD(1,ITERM) NASR FINVCG=FREAD(1,ITERM) NASR C NASR C COMPUTE THE OTHER PARAMETERS OF THE 'GEOCENTRIC' ELLIPSOID NASR C NASR FCG=1./FINVCG NASR BCG=ACG*(1.-FCG) NASR ECG=DSQRT(2.*FCG-FCG*FCG) NASR 5094 CONTINUE NASR C NASR C PRINT THE PARAMETERS AND TRANSLATIONS OF THE GEODETIC REFERENCE NASR C ELLIPSOID. NASR C NASR PRINT5080 NASR 5080 FORMAT('1') NASR PRINT592 NASR 592 FORMAT(///) NASR PRINT5003 NASR C C 5003 FORMAT(10X,'PARAMETERS OF THE GEODETIC REFERENCE ELLIPSOID :'/10X,NASR 1'----------------------------------------------') NASR PRINT30,GELLIP,A,B,FINV,E NASR 30 FORMAT(62X,'NAME = ',A3/65X,'A = ',F10.2,3X,'METERS'/65X,'B = ',F1NASR 10.2,3X,'METERS'/65X,'F = 1/',F10.6/65X,'E = ',F10.8) NASR PRINT592 NASR PRINT40,X0,Y0,Z0 NASR 40 FORMAT(10X,'GEODETIC DATUM TRANSLATIONS FROM THE GEOCENTRE :'/10X,NASR 1'----------------------------------------------'/60X,'XTRANS = ',FNASR 210.2,3X,'METERS'/60X,'YTRANS = ',F10.2,3X,'METERS'/60X,'ZTRANS = 'NASR 3,F10.2,3X,'METERS') NASR IF(CGELIP.EQ.ELNAME) GO TO 5084 NASR C NASR C PRINT THE PARAMETERS OF THE ARBITRARILY CHOZEN GEOCENTRIC ELLIPSOID NASR C NASR PRINT592 NASR PRINT5004 NASR 5004 FORMAT(10X,'PARAMETERS OF THE CHOZEN GEOCENTRIC ELLIPSOID :'/10X,'NASR 1---------------------------------------------') NASR PRINT30,CGELIP,ACG,BCG,FINVCG,ECG NASR 5084 CONTINUE NASR C NASR C***********************************************************************NASR C NASR C READ-IN THE CHOSEN STATION NUMBER (ICHOOZ) ,THIS STATION WILL BE THE NASR C ONLY STATION FOR WHICH THE PROGRAM WILL PRINT-OUT THE SUMMARY ,AND NASR C PUNCH (IF DESIRED) THE DATA NECESSARY TO PERFORM THE L.S. SOLUTION NASR C FOR THE COORDINATES OF THIS STATION (ORBIT-KNOWN CASE) NASR C NASR ICHOOZ=FREAD(1,ITERM) NASR C NASR C READ-IN THE CHOZEN TRACKING STATION ABBREVIATED NAME (STNAME) ,WHICH NASR C SHOULD BE 4 ALPHAMERIC CHARACTERS STARTING BY ALPHABETIC ; AND THE NASR C DIRECTION OF RECKONING THE LONGITUDE (LONDIR) ,WHICH SHOULD BE NASR C EITHER ; EAST OR WEST NAS9 C NASR READ5000,STNAME,LONDIR NASR 5000 FORMAT(A4,2X,A4) NASR C NASR C READ-IN THE APPROXIMATE GEOGRAPHICAL COORDINATES OF THE CHOZEN TRACK-NASR C ING STATION ; LATITUDE & LONGITUDE IN (DEG,MIN & SEC) ,AND THE GEOME-NASR C TRICAL HEIGHT ABOVE THE GEODETIC ELLIPSOID (IN METERS). NASR C NOTE THAT ALL PUNCHED VALUES FOR THE LONGITUDE SHOULD BE +VE NASR C NASR DLAT=FREAD(1,ITERM) NASR XLAT=FREAD(1,ITERM) NASR SLAT=FREAD(1,ITERM) NASR DLON=FREAD(1,ITERM) NASR XLON=FREAD(1,ITERM) NASR SLON=FREAD(1,ITERM) NASR OBSHT=FREAD(1,ITERM) NASR C NASR C COVERT THE APPROXIMATE (LAT ,LONG) TO RADIANS NASR C NASR OBSLAT=RADIAN(DLAT,XLAT,SLAT) NASR OBSLON=RADIAN(DLON,XLON,SLON) NASR IF(LONDIR.EQ.NEGTIV) OBSLON=-1.*OBSLON NASR C C C PRINT THE APPROX. COORDINATES (LAT,LONG)-IN DEG,MIN & SEC - AND THE NASR C HEIGHT IN METERS NASR C NASR PRINT592 NASR PRINT 41 NASR 41 FORMAT(' '/15X,'TABLE OF APPROXIMATE GEOGRAPHICAL COORDINATES OF NASR 1THE TRACKING STATION'/15X,'=======================================NASR 2=============================='/) NASR PRINT42 NASR 42 FORMAT(3X,'STATION',2X,'STATION',5X,'LATITUDE',10X,'LONGITUDE',11XNASR 1,'HEIGHT',4X,' ',5X,'DESCRIPTION'/5X,'NAME',4X,'NUMBER',3NASR 2X,'(+VE NORTH)',9X,'(+VE EAST)',9X,'(METERS)',3X,' '/21X,NASR 3'DEG',2X,'MIN',2X,'SEC',6X,'DEG',2X,'MIN',2X,'SEC') NASR IDEG1=IDEG(OBSLAT,MINT,SEC1) NASR IDEG2=IDEG(OBSLON,MIN2,SEC2) NASR PRINT43,STNAME,ICHOOZ,IDEG1,MINT,SEC1,IDEG2,MIN2,SEC2,OBSHT NASR 43 FORMAT(5X,A4,5X,I2,5X,I3,2X,I2,2X,F7.4,3X,I3,2X,I2,2X,F7.4,F13.5) NASR PRINT44,GELLIP NASR 44 FORMAT('+',82X,'*** GEODETIC-',A3,' ***') NASR C NASR C CALLING THE SUBROUTINE 'GEOREC' TO CONVERT THE APPROXIMATE GEOGRAPH- NASR C ICAL COORDINATES OF THE TRACKING STATION TO ITS CORRESPONDING GEO- NASR C CENTRIC (X,Y,Z) CARTESIAN COORDINATES ,AND PRINT THE LATER. NASR C NASR CALL GEOREC(A,B,X0,Y0,Z0,OBSLAT,OBSLON,OBSHT,OBSX,OBSY,OBSZ) NASR C NASR PRINT592 NASR PRINT45 NASR 45 FORMAT(' '/17X,'TABLE OF APPROXIMATE CARTESIAN COORDINATES OF THENASR 1 TRACKING STATION'/17X,'==========================================NASR 2========================'/) NASR PRINT46 NASR 46 FORMAT(5X,'STATION NAME',2X,'ST. NO.',6X,'X -(METERS)',6X,'Y -(METNASR 1ERS)',6X,'Z -(METERS)',10X,'DESCRIPTION'/) NASR PRINT47,STNAME,ICHOOZ,OBSX,OBSY,OBSZ NASR 47 FORMAT(9X,A4,8X,I2,3X,3(5X,F12.3)) NASR PRINT49,GELLIP NASR 49 FORMAT('+',82X,'***GEOCENTRIC-',A3,' ***') NASR C NASR C***********************************************************************NASR C NASR C READ-IN THE YEAR OF OBSERVATIONS ,AND THE APRIORI VARIANCE-FACTOR NASR C (IN CYCLES/SEC SQUARED) NASR C NASR IYEAR=FREAD(1,ITERM) NASR IAVF=FREAD(1,ITERM) NASR C NASR C INITIALIZATION OF SOME VARIABLES NASR C NASR C KPASS : IS THE TOTAL NUMBER OF THE ACCEPTED PASSES NASR C NDOPT : IS THE ACCUMULATED NUMBER OF OBSERVED SATELLITE POSITIONS NASR C NACDOP : IS THE ACCUMULATED NUMBER OF ACCEPTED DOPPLER COUNTS NASR C NASR KPASS=0 NASR NDOPT=0 NASR NACDOP=0 NASR C C C***********************************************************************NASR C NASR C RETURN HERE TO START READING-IN AND PROCESSING THE DATA OF A NEW PASSNASR C NASR C NOTE THAT THE FOLLOWING DELIMITERS ARE THE 'END CODES' OF THE DIFFER-NASR C ENT 'INPUT DATA' ,WHICH WILL ORGANIZE THE SEQUENCE OF THE PROGRAM NASR C STEPS BEING EXECUTED : NASR C 999 : IS A CODE INDICATING THE END OF THE SETS OF VARIABLE ORBITAL NASR C PARAMETERS FOR EACH INDIVIDUAL PASS. NASR C 99 : IS A CODE INDICATING THE END OF ALL DATA PERTAINING TO EACH NASR C INDIVIDUAL PASS ;THIS MEANS THE END OF ALL DOPPLER DATA TAKEN NASR C FROM ALL OBSERVING STATIONS TO THIS PARTICULAR PASS. NASR C 9999 : IS A CODE INDICATING THE END OF ALL INPUTTED DATA FOR ALL NASR C PASSES BEING PROCESSED. NASR C NASR 1000 CONTINUE NASR C NASR C***********************************************************************NASR C NASR C READ-IN THE SATELLITE IDENTIFICATION CARD ,WHICH CONSISTS OF : DAY # NASR C ,PASS # (W.R.T. THE DAY #) ,SATELLITE # ,AND LOCKON-TIME (HR & MIN) NASR C NASR IDAY=FREAD(1,ITERM) NASR C NASR C CHECK FOR THE END OF THE DATA ; IF YES ,JUMP TO PRINT THE SUMMARY OF NASR C THE OUTPUT FROM THIS PROGRAM. NASR C NASR IF(IDAY.EQ.9999) GOTO 990 NASR C NASR IPASS=FREAD(1,ITERM) NASR ISAT=FREAD(1,ITERM) NASR NHR=FREAD(1,ITERM) NASR MIN=FREAD(1,ITERM) NASR PRINT1001,KSEQ 5 KSEQ=KPASS+1 NASR 1001 FORMAT('1'/20X,'SATELLITE MAJORITY VOTING INPUT DATA *** PASS SEQUNASR 1ENCE NO.',I4,' ***'/20X,'===================================='/) NASR PRINT 940,IDAY,IPASS,ISAT,NHR,MIN NASR 940 FORMAT(10X,4HDAY=,I5,5X,5HPASS=,I5,5X,4HSAT=,I5,5X,5HHOUR=,I5,5X, NASR 1 4HMIN=,I5) NASR MIN1=60*NHR+MIN+0.5 NASR KTCX=MIN1-4*(MIN1/4) NASR C NASR C***********************************************************************NASR C NASR C READ-IN AND PRINT-OUT THE UNSCALED FIXED (OR MEAN) ORBITAL PARAMETERSNASR C NASR PRINT2052 NASR PRINT910 NASR 910 FORMAT(/5X,'UNSCALED SATELLITE MEAN ORBIT PARAMETERS'/5X,'--------NASR 1--------------------------------'/) NASR TP=FREAD(2,ITERM)/1E4 NASR TP=TP+DTPERG NASR XN=DEGRA*(3.+FREAD(2,ITERM)/1E7) NASR WNOT=DEGRA*(FREAD(2,ITERM)/1E4) NASR WDOT=DEGRA*(FREAD(2,ITERM)/1E7) NASR C C ECC=FREAD(2,ITERM)/1E6 NASR AO=FREAD( 2,ITERM)*10. NASR RA=DEGRA*(FREAD(2,ITERM)/1E4) NASR RADOT=DEGRA*(FREAD(2,ITERM)/1E7) NASR CI=FREAD(2,ITERM)/1E6 NASR GLAM=DEGRA*(FREAD(2,ITERM)/1E4) NASR SI=FREAD(2,ITERM)/1E6 NASR C NASR C COMPUTE THE TIME DIFFERENCE BETWEEN THE FIRST TIME OF PERIGEE AFTER NASR C THE LAST DATA INJECTION ,AND THE LOCKON-TIME. NASR C NASR DTP=MIN1-TP NASR IF(DTP.LE.0) DTP=DTP+1440. NASR IF(DTP+2.*PI/XN.GE.1440.) DTP=DTP-1440. NASR C NASR C***********************************************************************NASR C NASR C NASR C READ-IN AND PRINT THE UNSCALED VARIABLE (OR EPHEMERAL) ORBITAL NASR C PARAMETERS ; UP TO 17 SETS CAN BE READ ,WHICH SHOULD BE TERMINATED NASR C BY A DELIMITER '999' . NASR C NASR PRINT 915 NASR 915 FORMAT(//10X,'UNSCALED CORRECTIONS TO THE MEAN ORBIT PARAMETERS'/1NASR 10X,'-------------------------------------------------'/) NASR PRINT 501 NASR 501 FORMAT(5X,'TIME',4X,'ECCENTRIC-ANOMALY',3X,'SEMI-MAJOR AXIS',3X,'ONASR 1UT-OF-PLANE'/3X,'(UT MIN)',5X,'CORRECTION',9X,'CORRECTION',7X,'COMNASR 2PONENT'/) NASR DO 930 I=1,17 NASR DEK(I)= FREAD( 1,ITERM) NASR IF(DEK(I).EQ.999.) GOTO 119 NASR DAK(I)= FREAD( 1,ITERM) NASR ETAK(I)= FREAD( 1,ITERM) NASR MINK(I)=MIN+2*(I-2)+0.5 NASR IF(MINK(I).LT.0) MINK(I)=MINK(I)+60 NASR IF(MINK(I).GE.60) MINK(I)=MINK(I)-60 NASR NUM=I-2 NASR PRINT 931,MINK(I),DEK(I),DAK(I),ETAK(I) NASR 931 FORMAT(4X,I4,11X,F5.0,14X,F5.0,11X,F5.0) NASR C NASR C SCALING THE CORRECTIONS TO THE MEAN ORBIT PARAMETERS NASR C NASR DEK(I)=DEGRA*(DEK(I)/1E4) NASR DAK(I)=DAK(I)*10. NASR 930 ETAK(I)=ETAK(I)*10. NASR C NASR C***********************************************************************NASR C NASR C PERFORMING THE LAGRANGIAN INTERPOLATION FOR THE OUT-OF-PLANE COMPONENTASR C NASR 119 IF(KTCX.EQ.0) GO TO 118 NASR NUM1=NUM+1 NASR DO 120 I=2,NUM1,2 NASR ETAK(I)=(ETAK(I-1)+ETAK(I+1))/2. NASR J=I-1 NASR IF(I.EQ.NUM1) J=I-3 NASR 120 ETAK(I)=ETAK(I)-(ETAK(J)+ETAK(J+4)-2.*ETAK(J+2))/8. NASR GO TO 122 NASR C C C 118 DO 121 I=3,NUM,2 NASR ETAK(I)=(ETAK(I-1)+ETAK(I+1))/2. NASR J=I-1 NASR IF(I.EQ.NUM) J=I-3 NASR 121 ETAK(I)=ETAK(I)-(ETAK(J)+ETAK(J+4)-2.*ETAK(J+2))/8. NASR 122 CONTINUE NASR C NASR C ELIMINATE THE FIRST AND LAST ROWS OF THE VARIABLE PARAMETERS TABLE NASR C ,TO REDUCE THEIR NUMBER TO THE NUMBER OF OBSERVED SATELLITE POSITIONSNASR C NASR DO 400 I=1,NUM NASR MINK(I)=MINK(I+1) NASR DEK(I)=DEK(I+1) NASR DAK(I)=DAK(I+1) NASR 400 ETAK(I)=ETAK(I+1) NASR C NASR C***********************************************************************NASR C NASR C PRINT-OUT THE SCALED MEAN ORBIT PARAMETERS ,AND THE SCALED VARIABLE NASR C PARAMETERS (CORRECTIONS TO SOME OF THE MEAN ORBIT PARAMETERS) NASR C NASR PRINT 502 NASR 502 FORMAT('1'//5X,'SCALED MEAN ORBIT PARAMETERS : ( UNITS ARE IN :- NASR 1METERS ,RADIANS ,AND UT MINUTES )'/5X,'---------------------------NASR 2---'/) NASR PRINT861,TP,XN NASR 861 FORMAT(5X,F16.10,2X,': IS THE TIME OF PERIGEE -(UT MINUTES)'//5X,FNASR 116.10,2X,': IS THE SATELLITE MEAN MOTION -(RADIANS/MIN)'/) NASR PRINT503,WNOT,WDOT NASR 503 FORMAT(5X,F16.10,2X,': IS THE ARGUMENT OF PERIGEE AT THE TIME OF PNASR 1ERIGEE -(RADIANS)'//5X,F16.10,2X,': IS THE PRECESSION RATE OF PERINASR 2GEE -(RADIANS/MIN)'/) NASR PRINT504,ECC,AO NASR 504 FORMAT(5X,F16.10,2X,': IS THE ECCENTRICITY OF THE ORBITAL ELLIPSE NASR 1-(UNITLESS)'//2X,F19.10,2X,': IS THE MEAN SEMI-MAJOR AXIS OF THE ONASR 2RBITAL ELLIPSE -(METERS)'/) NASR PRINT505,RA,RADOT NASR 505 FORMAT(5X,F16.10,2X,': IS THE RIGHT ASCENSION OF THE ASCENDING NODNASR 1E AT THE TIME OF PERIGEE -(RADIANS)'//5X,F16.10,2X,': IS THE PRECENASR 2SSION RATE OF THE ASCENDING NODE -(RADIANS/MIN)'/) NASR PRINT506,CI,GLAM,SI NASR 506 FORMAT(5X,F16.10,2X,': IS THE COSINE OF THE ORBIT INCLINATION -(UNNASR 1ITLESS)'//5X,F16.10,2X,': IS THE GREENWICH APPARENT SIDERIAL TIME NASR 2AT THE TIME OF PERIGEE -(RADIANS)'//5X,F16.10,2X,': IS THE SINE OFNASR 3 THE ORBIT INCLINATION -(UNITLESS)'/) NASR PRINT862 NASR 862 FORMAT(///20X,'SCALED VARIABLE PARAMETERS'/20X,'------------------NASR 1--------') NASR PRINT507 NASR 507 FORMAT(/,5X,'TIME',7X,'ECCENTRIC-ANOMALY',6X,'SEMI-MAJOR AXIS',6X,NASR 1'OUT-OF-PLANE'/3X,'(UT MIN)',3X,'CORRECTION (RADIANS)',3X,'CORRECTNASR 1ION (METERS)',2X,'COMPONENT (METERS)'/) NASR DO 865 I=1,NUM NASR 865 PRINT 863,MINK(I),DEK(I),DAK(I),ETAK(I) NASR 863 FORMAT(4X,I4,10X,F13.10,10X,F11.4,9X,F9.4/) NASR C NASR C***********************************************************************NASR C NASR C C C COMPUTE THE SATELLITE (X,Y & Z) CARTESIAN COORDINATES ,REFERRED TO NASR C THE AVERAGE TERRESTRIAL COORDINATE SYSTEM. NASR C NASR C INPUTS NASR C DTP=TIME BETWEEN TIME OF PERIGEE AND LOCKON TIME (MINUTES) NASR C K = NUMBER OF POSITIONS TO BE COMPUTED NASR C TNK = TIME INTERVAL FOR DOPPLER COUNT NK=2 MIN IN DTK EQN BELOWNASR C ALL SATELLITE ORBIT PARAMETERS FROM MESSAGE NASR C OUTPUT NASR C SATELLITE COORDINATES AT TWO MINUTE INTERVALS (XS,YS,ZS) NASR C REFERRED TO THE AVERAGE TERRESTRIAL COORDINATE SYSTE. NASR C NASR C COMPUTE THE SATELLITE (XYZ) CARTESIAN COORDINATES ,W.R.T. ITS ORBIT NASR C COORDINATE SYSTEM NASR C NASR PRINT 210 NASR 210 FORMAT('1'//10X,'TABLE OF RECTANGULAR COORDINATES OF SATELLITE POSNASR 1ITIONS REFERRED TO ITS ORBITAL SYSTEM'/10X,'======================NASR 1================================================================')NASR PRINT508 NASR 508 FORMAT(/5X,'POSITION',5X,'TIME',14X,'X',18X,'Y',18X,'Z',11X,'GEOCENASR 1NTRIC RANGE'/6X,'NUMBER',4X,'(UT MIN)',8X,'(METERS)',11X,'(METERS)NASR 2',11X,'(METERS)',11X,'(METERS)'/) NASR DO 250 K=1,NUM NASR DTK=DTP+2.*(K-1.) NASR FMK=XN*DTK NASR EK=FMK+ECC*DSIN(FMK)+DEK(K) NASR AK=AO+DAK(K) NASR X(1)=AK*(DCOS(EK)-ECC) NASR X(2)=AK*DSIN(EK) NASR X(3)=ETAK(K) NASR X(4)=DSQRT(X(1)*X(1)+X(2)*X(2)+X(3)*X(3)) NASR PRINT220,K,MINK(K),(X(I),I=1,4) NASR 220 FORMAT(8X,I1,10X,I2,3X,4(6X,F13.4)) NASR C NASR C TRANSFORMATION OF THE SATELLITE (XYZ) COORDINATES FROM THE ORBITAL NASR C SYSTEM TO THE AVERAGE TERRESTRIAL SYSTEM (VIA THE APPARENT CELESTIAL NASR C SYSTEM) ; USING THE ORTHOGONAL ROTATION MATRICES. NASR C NASR DATA NAXIS, NROT/3,1,3,3/ NASR ANGLE(1)=-(WNOT-WDOT*DTK) NASR ANGLE(2)=-DATAN2(SI,CI) NASR ANGLE(3)=-(RA+RADOT*DTK-(GLAM+WE*DTK)) NASR C NASR C CALLING THE SUBROUTINE 'DROTAT' TO EVALUATE THE PRODUCT ROTATION NASR C MATRIX (ROT) ,NECESSARY TO PERFORM THE COORDINATE TRANSFORMATION NASR C NASR CALL DROTAT(ANGLE,NAXIS,NROT,ROT) NASR C NASR C COMPUTE THE SATELLITE (XYZ) AVERAGE TERRESTRIAL ,BY MULTIPLYING THE NASR C CORRESPONDING ORBITAL (XYZ) BY THE PRODUCT ROTATION MATRIX (ROT) NASR C XYZ(TERR)=ROT*XYZ(ORBIT) NASR C NASR XS(K,1)=0.0 NASR XS(K,2)=0.0 NASR XS(K,3)=0.0 NASR C C DO 230 I=1,3 NASR DO 230 J=1,3 NASR 230 XS(K,I)=XS(K,I)+X(J)*ROT(I,J) NASR XS(K,4)=DSQRT(XS(K,1)*XS(K,1)+XS(K,2)*XS(K,2)+XS(K,3)*XS(K,3)) NASR C NASR C***********************************************************************NASR C NASR C THE SATELLITE POSITIONS REFERRED TO THE TERRESTRIAL REFERENCE ELLIPS-NASR C OID. NOTE THAT THE COMPUTED (LAT,LONG) IN THIS CASE WILL DETERMINE NASR C COMPUTE THE GEODETIC COORDINATES (LATITUDE , LONGITUDE & HEIGHT) OF NASR C THE POSITION OF THE SATELLITE SUBPOINTS ,WHICH ARE THE TRACING OF NASR C THE SATELLITE TRACK ON THE EARTH ELLIPSOID SURFACE (BY USING THE NASR C ORTHOGONAL PROJECTION) NASR C NASR SATX=XS(K,1) NASR SATY=XS(K,2) NASR SATZ=XS(K,3) NASR C NASR C CALLING THE SUBROUTINE 'RECGEO' TO CONVERT THE (XYZ) CARTESIAN COORD-NASR C INATES OF THE SATELLITE POSITIONS TO ITS CORRESPONDING (LAT,LONG,HT) NASR C W.R.T. THE EARTH REFERRENCE ELLIPSOID. NASR C NASR CALL RECGEO(A,B,X0,Y0,Z0,SATX,SATY,SATZ,SUBPHI,SUBLAM,SUBHT) NASR SUB(K,1)=SUBPHI NASR SUB(K,2)=SUBLAM NASR SUB(K,3)=SUBHT NASR C NASR 250 CONTINUE NASR C NASR C***********************************************************************NASR C NASR C PRINT THE SATELLITE (XYZ) AVERAGE TERRESTRIAL CARTESIAN COORDINATES NASR C NASR PRINT259 NASR 259 FORMAT(//10X,'TABLE OF SATELLITE RECTANGULAR COORDINATES REFEERED NASR 1TO THE AVERAGE TERRESTRIAL SYSTEM'/10X,'==========================NASR 2==========================================================='/) NASR PRINT508 NASR DO 260 K=1,NUM NASR PRINT220,K,MINK(K),(XS(K,I),I=1,4) NASR 260 CONTINUE NASR C NASR C***********************************************************************NASR C NASR C PRINT THE SATELLITE SUBPOINT COORDINATES (LAT,LONG) AND THE SATELLITENASR C HEIGHT W.R.T. THE EARTH REFERENCE ELLIPSOID. NASR C NASR PRINT 510 NASR 510 FORMAT(' '//,15X,'TABLE OF COMPUTED GEODETIC COORDINATES OF THE SANASR 1TELLITE POSITIONS'/15X,'==========================================NASR 2======================='/) NASR PRINT 511 NASR 511 FORMAT(10X,'POSITION',5X,'TIME',8X,'LATITUDE',11X,'LONGITUDE',13X,NASR 1'HEIGHT'/11X,'NUMBER',4X,'(UT MIN)',4X,'(+VE NORTH)',10X,'(+VE EASNASR 2T)',11X,'(METERS)'/32X,'DEG',2X,'MIN',2X,'SEC',7X,'DEG',2X,'MIN',2NASR 3X,'SEC') NASR C C DO 270 K=1,NUM NASR RLAT=SUB(K,1) NASR RLON=SUB(K,2) NASR LATD=IDEG(RLAT,LATM,SLAT) NASR LOND=IDEG(RLON,LONM,SLON) NASR 270 PRINT265,K,MINK(K),LATD,LATM,SLAT,LOND,LONM,SLON,SUB(K,3) NASR 265 FORMAT(13X,I1,10X,I2,1X,2(4X,I3,2X,I2,2X,F7.4),5X,F12.4) NASR C NASR C***********************************************************************NASR C NASR C READ-IN THE TRACKING STATION-NUMBER ,FOLLOWED BY THE SET OF DOPPLERS NASR C OBSERVED FROM THAT STATION. NOTE THAT **ANY** STATION-NUMBERS WITH NASR C THEIR DOPPLER DATA FROM THE SAME PASS CAN BE READ-IN HERE ,BUT EVEN- NASR C TUALLY ALL OF THEM WILL BE IGNORED EXCEPT THE 'CHOZEN STATION-ICHOOZ'NASR C THIS IS MOTIVATED BY THE POSSIBILITY OF USING THE SAME DATA AVAILABLENASR C FOR THE SIMULTANEOUS FORMAT FROM MANY STATIONS ,WITHOUT REARRANGE- NASR C MENT SPECIFICALLY FOR THIS PROGRAM. AGAIN THEY HAVE NO EFFECT HERE, NASR C NASR C THE REASON FOR USING TWO DIFFERENT NAMES FOR THE STATION-NUMBER IS NASR C TO BE SURE AFTER PRINTING THE RESULTS THAT THE DOPPLER DATA BEING NASR C PROCESSED PERTAINS TO THE CHOZEN TRACKING STATION 'ICHOOZ' ,AND NOT NASR C FOR ANY OTHER OBSERVING STATION. NASR C NASR 900 ISTN=FREAD(1,ITERM) NASR C NASR C CHECK FOR THE END OF THE PASS (CODE '99') ,IF YES ;JUMP TO READ AND NASR C PROCESS DATA FOR A NEW PASS (IF ANY) NASR C NASR IF(ISTN.EQ.99) GOTO 901 NASR IF(ISTN.EQ.ICHOOZ) GOTO 340 NASR C NASR C***********************************************************************NASR C NASR C READ-IN AND IGNORE ANY EXISTING DOPPLER DATA INPUTTED FROM ANY STAT- NASR C ION ,OTHER THAN THE CHOZEN STATION 'ICHOOZ' NASR C NASR DO 330 KK=2,NUM NASR DUMMY=FREAD(1,ITERM) NASR 330 DUMMY=FREAD(1,ITERM) NASR GOTO 900 NASR C NASR C***********************************************************************NASR C NASR 340 CONTINUE NASR PRINT351,STNAME,ISTN NASR 351 FORMAT('1'///15X,'THE FOLLOWING RESULTS WERE COMPUTED FROM ** ',A4NASR 1,'-STATION ** NUMBER : ',I2) NASR C NASR C***********************************************************************NASR C NASR C COMPUTE THE TOPOCENTRIC RANGE VECTOR BETWEEN THE TRACKING STATION NASR C AND SATELLITE POSITIONS REFERRED TO THE AVERAGE TERRESTRIAL SYSTEM. NASR C NASR DO 610 K=1,NUM NASR R(K,1)=XS(K,1)-OBSX NASR R(K,2)=XS(K,2)-OBSY NASR C C R(K,3)=XS(K,3)-OBSZ NASR R(K,4)=DSQRT(R(K,1)*R(K,1)+R(K,2)*R(K,2)+R(K,3)*R(K,3)) NASR C NASR C TRANSFORMATION OF THE TOPOCENTRIC RANGE VECTOR FROM THE AVERAGE TERR.NASR C TO THE TOPOCENTRIC GEODETIC COORDINATE SYSTEM ; TO BE USED IN COMPUT-NASR C ING THE TOPOCENTRIC DIRECTIONS (AZIMUTH & ELEVATION ANGLE) FROM THE NASR C TRACKING STATION TO EACH SATELLITE POSITION. NASR C NASR DATA NAXTOP,NUMTOP/3,2,2/ NASR TOPANG(1)=OBSLON-PI NASR TOPANG(2)=OBSLAT-PI/2. NASR C NASR C CALLING THE SUBROUTINE 'DROTAT' TO EVALUATE THE REQUIRED PRODUCT NASR C ROTATION MATRIX (TOPROT) NEEDED FOR THE ABOVE TRANSFORMATION : NASR C 'TOPOCENTRIC-VECTOR' = TOPROT * 'AVER. TERR. -VECTOR' NASR C NASR CALL DROTAT (TOPANG,NAXTOP,NUMTOP,TOPROT) NASR C NASR RT(K,1)=0. NASR RT(K,2)=0. NASR RT(K,3)=0. NASR DO 235 I=1,3 NASR DO 235 J=1,3 NASR 235 RT(K,I)=RT(K,I)+R(K,J)*TOPROT(I,J) NASR RT(K,4)=DSQRT(RT(K,1)*RT(K,1)+RT(K,2)*RT(K,2)+RT(K,3)*RT(K,3)) NASR C NASR C REFLECT THE SECONDARY AXIS (RT(K,2)) TO OBTAIN A LEFT-HANDED SYSTEM NASR C ,WHICH IS ONE OF THE CHARACTERISTICS OF THE TOPOCENTRIC COORD. SYSTEMNASR C NASR RT(K,2)=-RT(K,2) NASR C NASR C COMPUTE THE AZIMUTH OF SATELLITE POSITIONS AS VIEWED FROM OBSERVER. NASR C NASR AZ(K)=DATAN2(RT(K,2),RT(K,1)) NASR C NASR C COMPUTE THE ELEVATION ANGLE OF EACH SATELLITE POSITION W.R.T. THE NASR C TRACKING STATION. NASR C NASR P=DSQRT(RT(K,1)*RT(K,1)+RT(K,2)*RT(K,2)) NASR VA(K)=DATAN2(RT(K,3),P) NASR 610 CONTINUE NASR C NASR C***********************************************************************NASR C NASR C PRINT-OUT THE SATELLITE RECTANGULAR COORDINATES REFERRED TO THE NASR C LOCAL GEODETIC COORDINATE SYSTEM. NASR C NASR PRINT 540 NASR 540 FORMAT(////10X,'TABLE OF SATELLITE RECTANGULAR COORDINATES REFERRENASR 1D TO THE LOCAL GEODETIC SYSTEM'/10X,'=============================NASR 2==================================================='/) NASR PRINT 541 NASR 541 FORMAT(5X,'POSITION',5X,'TIME',14X,'X',18X,'Y',18X,'Z',10X,'TOPOCENASR 1NTRIC RANGE'/6X,'NUMBER',4X,'(UT MIN)',8X,'(METERS)',11X,'(METERS)NASR 2',11X,'(METERS)',11X,'(METERS)'/) NASR DO 280 K=1,NUM NASR C C 280 PRINT220,K,MINK(K),(RT(K,I),I=1,4) NASR C NASR C PRINT THE TOPOCENTRIC RANGE AND TOPOCENTRIC DIRECTIONS (AZIMUTH & NASR C ELEVATION ANGLE) FROM THE TRACKING STATION TO THE SATELLITE POSITIONSNASR C NASR PRINT592 NASR PRINT 291 NASR 291 FORMAT(' '//10X,'TABLE OF COMPUTED ; TOPOCENTRIC RANGES , AZIMUTH NASR 1AND VERTICAL ANGLE TO THE SATELLITE'/10X,'========================NASR 2============================================================'/) NASR PRINT 543 NASR 543 FORMAT(10X,'POSITION',6X,'TIME',5X,'TOPOCENTRIC-RANGE',5X,'AZIMUTHNASR 1',10X,'VERTICAL-ANGLE'/11X,'NUMBER',5X,'(UT MIN)',7X,'(METERS)',8XNASR 2,'DEG',2X,'MIN',2X,'SEC',7X,'DEG',2X,'MIN',2X,'SEC'/) NASR DO 290 K=1,NUM NASR PAZ=AZ(K) NASR PVA=VA(K) NASR IDAZ=IDEG(PAZ,MAZ,SAZ) NASR IDVA=IDEG(PVA,MVA,SVA) NASR 290 PRINT285,K,MINK(K),R(K,4),IDAZ,MAZ,SAZ,IDVA,MVA,SVA NASR 285 FORMAT(13X,I2,10X,I2,8X,F13.4,2(4X,I4,2X,I2,2X,F6.3)) NASR C NASR C PRINT THE SUMMARY OF IMPORTANT INFORMATIONS FOR THE SATELLITE PASS NASR C AT THE TIME OF ITS CLOSEST APPROACH ,NAMELY ; THE MAX. ELEVATION ,THENASR C DIRECTION OF THE TRACK ,AND THE AZIMUTH - W.R.T. THE TRACKING STATIONNASR C NASR DATA INOR,ISOU,IEAST,IWEST/'S-N','N-S','EAST','WEST'/ NASR DO 611 K=2,NUM NASR IF(VA(K).LT.VA(K-1)) GOTO 612 NASR 611 CONTINUE NASR 612 VAMAX=VA(K-1)-(VA(K)-VA(K-2))**2/(8.*(VA(K)-2.*VA(K-1)+VA(K-2))) NASR VAMAX=VAMAX*180./PI NASR IDIR=INOR NASR IAZ=IEAST NASR IF(SUB(K,1).LT.SUB(K-1,1)) IDIR=ISOU NASR IF(AZ(K).LT.0.) IAZ=IWEST NASR PRINT613,VAMAX,IDIR,IAZ NASR 613 FORMAT('1'//5X,'SUMMARY OF INFORMATIONS AT THE SATELLITE CLOSEST ANASR 1PPROACH :'/5X,'---------------------------------------------------NASR 2-----'/65X,'MAX. ELEVATION =',F7.3,3X,'DEGREES'/65X,'TRACK DIRECTNASR 3ION =',A6/65X,'TRACK AZIMUTH =',A6/) NASR C NASR C***********************************************************************NASR C NASR C READ-IN ,ANALIZE AND PRINT THE DOPPLER DATA ; EITHER 'MAGNAVOX' OR NASR C 'ITT' FORMAT ,WOULD BE ACCEPTED. NASR C NASR PRINT 631 NASR 631 FORMAT(//,25X,'DETAILS OF THE OBSERVED AND CORRECTED DOPPLER DATA NASR 1; UNITS ARE IN CYCLES'/25X,'======================================NASR 2=================================='/) NASR PRINT 545 NASR 545 FORMAT(3X,'TIME',3X,'ELEVATION-ANGLE',3X,'DOPPLERS',7X,'OBSERVED-DNASR 1OPPLERS',3X,'IONOSPHERIC-',3X,'TROPOSPHERIC',3X,'REDUCED-DOPPLERS'NASR 2/1X,'(UT MIN)',4X,'(DEGREES)',4X,'SEQUENCE NO.',4X,'400-MHZ',5X,'1NASR 350-MHZ',3X,'CORRECTION',5X,'CORRECTION',5X,'(TO THE VACIUM)'/) NASR C C PVA=VA(1)*180./PI NASR PRINT554,MINK(1),PVA NASR 554 FORMAT(4X,I2,6X,F8.3) NASR C NASR C INITIALIZE ZERO VALUES FOR THE DOPPLER DATA AND THIER CORRECTIONS. NASR C NASR IDOP(1)=0. NASR DOPC(1)=0. NASR REFR(1)=0. NASR DOPR(1)=0. NASR VDOP(1)=0. NASR DO 920 K=2,NUM NASR C NASR C READ-IN THE PAIRS OF DOPPLER COUNTS ,NAMELY ; THE '400 MHZ' AND THE NASR C '150 MHZ'-(REFRACTION) DOPPLERS. NASR C NASR DOPC(K)=FREAD(1,ITERM) NASR REFR(K)=FREAD(1,ITERM) NASR C NASR C IDENTIFY THE ZERO DOPPLERS ; IF EITHER THE 400-MHZ OR THE 150-MHZ NASR C DOPPLER EQUALS ZERO ,THEN ASSIGN ZERO VALUES FOR IONO. & TROPO. CRRE-NASR C CTIONS AND FOR THE FINAL REDUCED DOPPLER TO THE VACIUM AS WELL. NASR C NASR IDOP(K)=0 NASR DOPR(K)=0.D 00 NASR VDOP(K)=0.D 00 NASR DIFF=0.D0 NASR DIFF2=0.D0 NASR PVA=VA(K)*180./PI NASR IF(DOPC(K).EQ.0.OR.REFR(K).EQ.0) GO TO 632 NASR C NASR IDOP(K)=K-1 NASR C CORRECT THE OBSERVED DOPPLERS FOR THE IONOSPHERIC REFRACTION ,WHICH NASR C IS A FREQUENCY DEPENDENT ; WE CAN GET HERE A FIRST ORDER CORRECTION NASR C NASR C 1- FOR THE 'ITT'-FORMAT NASR C NASR IF(REFR(K).LT.5000.) DOPR(K)=DOPC(K)-(REFR(K)-2000.)*24./55. NASR C NASR C 2- FOR THE 'MAGNAVOX'-FORMAT NASR C NASR IF(REFR(K).GT.2000000.) DOPR(K)=DOPC(K)-(REFR(K)-DOPC(K))*9./55. NASR C NASR C APPLY THE TROPOSPHERIC REFRACTION CORRECTION TO THE CORRECTED DOPPLERNASR C DUE TO IONOSPHERE ,TO OBTAIN THE REDUCED DOPPLERS TO THE VACIUM, NASR C *** USING 'HOPFIELD'-MODEL 1969 *** NASR C NASR ELEVD1=DSQRT(VA(K-1)**2+THD) NASR ELEVD2=DSQRT(VA(K )**2+THD) NASR ELEVW1=DSQRT(VA(K-1)**2+THW) NASR ELEVW2=DSQRT(VA(K )**2+THW) NASR TROP2=XKD/DSIN(ELEVD2)+XKW/DSIN(ELEVW2) NASR TROP1=XKD/DSIN(ELEVD1)+XKW/DSIN(ELEVW1) NASR VDOP(K)=DOPR(K)-WAVENO*(TROP2-TROP1) NASR C NASR C PRINT-OUT THE TABLE OF OBSERVED AND CORRECTED DOPPLER DATA. NASR C C DIFF=DOPR(K)-DOPC(K) NASR DIFF2=VDOP(K)-DOPR(K) NASR 632 PRINT630,IDOP(K),DOPC(K),REFR(K),DIFF,DIFF2,VDOP(K),MINK(K),PVA NASR 630 FORMAT(30X,I2,8X,F10.1,2X,F10.1,2X,F8.3,7X,F8.3,8X,F12.3/4X,I2,6X,NASR 1F8.3) NASR 920 CONTINUE NASR PRINT2052 NASR C NASR C***********************************************************************NASR C NASR C REJECT HORIZON DOPPLERS AND COUNTING THE REMAINING USEFUL ONES. NAK9 C HOR : IS THE OPTIONAL HORIZON CRITERION IN DEGREES OF ELEVATION NASR C UNDER WHICH WE REJECT THE HORIZON DOPPLERS NASR C NREJEC : IS THE REJECTED NUMBER OF DOPPLERS UNDER THE ABOVE CRITERIONNASR C NDOP : IS THE NUMBER OF ACCEPTED DOPPLERS FOR EACH PASS. NASR C NASR HOR=HORREJ NASR NREJEC=0 NASR NDOP=0 NASR DO 710 K=1,NUM NASR IF(IDOP(K).EQ.0) GO TO 710 NASR IF(VA(K-1)*180./PI.GT.HOR) GO TO 703 NASR GO TO 704 NASR 703 IF(VA(K)*180./PI.GT.HOR) GO TO 705 NASR 704 IDOP(K)=0 NASR NREJEC=NREJEC+1 NASR GO TO 710 NASR 705 NDOP=NDOP+1 NASR 710 CONTINUE NASR PRINT7101,NREJEC,HOR NASR 7101 FORMAT('0',10X,I10,10X,'DOPPLERS REJECTED ON ',F5.1,5X,'DEGREES HONASR 1RIZON CRITERION') NASR C NASR C REJECT PASSES HAVING LESS THAN TWO-DOPPLERS (I.E. NDOP < 2 ) NASR C NASR PRINT2052 NASR PRINT712,IDAY,IPASS,ISAT NASR 712 FORMAT('0',3X,'PASS ON DAY',I5,5X,'PASS NUMBER',I5,5X,'SATELLITE NNASR 1UMBER',I5) NASR IF(NDOP.GE.2) GO TO 716 NASR PRINT714 NASR 714 FORMAT('+',68X,'**REJECTED-HAS LESS THAN TWO-DOPPLERS**') NASR GO TO 900 NASR C NASR C ACCEPT PASSES HAVING MORE THAN TWO-DOPPLERS (I.E. NDOP > 2 ) NASR C NASR 716 PRINT718 NASR 718 FORMAT('+',68X,'**ACCEPTED-HAS MORE THAN TWO-DOPPLERS**') NASR PRINT2052 NASR PRINT708,STNAME,ISTN,VAMAX,(IDOP(K),K=2,NUM),NDOP NASR 708 FORMAT('0',4X,'ST. NAME',3X,'ST. NO.',5X,'MAX. ELEV.',8X,'DOPPLERSNASR 1 TO BE USED',10X,'AVAILABLE DOPPLER EQUATIONS'//7X,A4,7X,I2,10X,F6NASR 2.3,7X,8I3,19X,I2) NASR C NASR C***********************************************************************NASR C NASR C C C STORE THE IMPORTANT INFORMATIONS ABOUT THE ACCEPTED PASSES ,WHICH NASR C WILL BE NEEDED ;EITHER IN THE L.S. SOLUTION TO DETERMINE THE COORD. NASR C OF THE TRACKING STATION ,OR THEY WILL BE RESERVED FOR FUTURE USE. NASR C NASR C STORE THE IDENTIFICATION CARD OF EACH ACCEPTED PASS. NASR C IDLSPE(KPASS) : IS THE I.D. CARD OF PASS # KPASS ,IT IS AN INTEGER NASR C NUMBER COMPOSED OF 12-DIGITS AND IT READS FROM LEFT NASR C TO RIGHT AS FOLLOWS : 3-DIGITS FOR THE DAY # ,2-DIGITSASR C FOR THE LOCKON HOUR ,2-DIGITS FOR THE LOCKON MINUTE NASR C ,2-DIGITS FOR THE SAT. # ,2-DIGITS FOR THE PASS # (W.NASR C R.T. THE DAY #) ,AND THE LAST DIGIT FOR THE ACCEPTED NASR C NUMBER OF DOPPLERS FROM THIS PASS # KPASS. NASR C NOTE THAT IDLSPE(12-DIGITS) IS A COMBINATION OF TWO ADJACENT INTEGERSNASR C VARIABLES :- IDL(7-DIGITS) & ISPE(5-DIGITS) ,THE 1ST GIVES ;DAY#/LOCKNASR C HR & MIN ,AND THE 2ND GIVES ;SAT#/PASS#/# OF DOPPLERS RESPECTIVELY , NASR C THE REASON FOR THIS SPLITTING IS THAT THE IBM/360 CANNOT HANDLE MORE NASR C THAN 10-DIGITS INTEGER CONSTANTS NASR C NASR KPASS=KPASS+1 NASR IDL(KPASS)=IDAY*10000+NHR*100+MIN NASR ISPE(KPASS)=ISAT*1000+IPASS*10+NDOP NASR C NASR C STORE THE CHARACTERISTICS OF EACH ACCEPTED PASS W.R.T. THE TRACKING NASR C STATION ,NAMELY ; THE MAX. ELEVATION ,THE DIRECTION OF THE TRACK ,THENASR C AZIMUTH ,THE QUADRANT ,AND THE NUMBER OF OBSERVED POSITIONS. NASR C NASR VAMSUM(KPASS)=VAMAX NASR IDRSUM(KPASS)=IDIR NASR IAZSUM(KPASS)=IAZ NASR DATA IZQUAD/'N-E','N-W','S-W','S-E'/ NASR IF(IAZ.EQ.IEAST.AND.IDIR.EQ.INOR) IQUADR(KPASS)=IZQUAD(1) NASR IF(IAZ.EQ.IEAST.AND.IDIR.EQ.ISOU) IQUADR(KPASS)=IZQUAD(4) NASR IF(IAZ.EQ.IWEST.AND.IDIR.EQ.INOR) IQUADR(KPASS)=IZQUAD(2) NASR IF(IAZ.EQ.IWEST.AND.IDIR.EQ.ISOU) IQUADR(KPASS)=IZQUAD(3) NASR NUMXYZ(KPASS)=NUM NASR C NASR C COMPUTE THE ROW-DIMENSION OF THE DESIGN MATRIX-'A' ,WHICH IS EQUAL NASR C TO THE ACCUMULATED NUMBER OF ACCEPTED DOPPLER COUNTS (I.E. EQUATIONS)NASR C NASR NACDOP=NACDOP+NDOP NASR NAROWS=NACDOP NASR C NASR C COMPUTE THE COLUMN-DIMENSION OF THE DESIGN MATRIX-'A' ,WHICH IS EQUALNASR C TO THE ACCUMULATED NUMBER OF THE UNKNOWN PARAMETERS (3 FOR THE TRACK-NASR C ING STATION COORD. + ONE FREQUENCY OFFSET PER EACH ACCEPTED PASS) NASR C NASR NACOLS=3+KPASS NASR C NASR C STORE THE SATELLITE COORDINATES AND CORRESPONDING DOPPLER COUNTS. NASR C (WHICH WERE TAKEN FROM THE CHOSEN TRACKING STATION ; ICHOOZ ) NASR C NASR DO 730 K=1,NUM NASR NDOPT=NDOPT+1 NASR IEQN(NDOPT)=IDOP(K) NASR DOP(NDOPT)=VDOP(K) NASR IF(IEQN(NDOPT).EQ.0) DOP(NDOPT)=0. NASR C C DO 730 J=1,3 NASR XSAT(NDOPT,J)=XS(K,J) NASR 730 CONTINUE NASR C NASR C RETURN-BACK TO READ AND IGNORE ALL SETS OF DOPPLER DATA FOR THE NASR C SAME PASS FROM ANY OTHER OBSERVING STATION (IF ANY) NASR C NASR GO TO 900 NASR C NASR C***********************************************************************NASR C NASR C RETURN-BACK TO READ AND PROCESS THE INPUT DATA FOR A NEW PASS. NASR C NASR 901 GOTO 1000 NASR C NASR C JUMP HERE TO PRINT A MESSAGE INDICATING THE END OF DATA FROM ALL NASR C PASSES ,AND TO START PRINTING THE REQUIRED OUTPUT-SUMMARY. NASR C NASR 990 CONTINUE NASR PRINT5080 NASR DO581 I=1,4 NASR 581 PRINT592 NASR PRINT993,KPASS NASR 993 FORMAT(/,18X,'**** E N D O F D A T A *** N O M O R E P A S NASR 1S E S ****'//38X,'========================='/,//38X,I3,2X,'PASSESNASR 2 WERE PROCESSED'//,/38X,'========================='/) NASR C NASR C***********************************************************************NASR C NASR C PRINT THE SUMMARY OF OUTPUT INFORMATIONS FROM THE 'PRE-BATCH' PROGRAMNASR C NOTE THAT THIS SUMMARY WILL BE PRINTED AS MANY AS (NCOPYS) NASR C NASR IREPET=0 NASR 2030 CONTINUE NASR PRINT2035 NASR 2035 FORMAT('1'//35X,'SUMMARY OF THE (PRE-BATCH)-PROGRAM OUTPUT'/35X,'=NASR 1========================================') NASR C NASR C NASR C PRINT THE PARAMETERS AND TRANSLATIONS OF THE GEODETIC ELLIPSOID NASR C NASR XTRANS=X0 NASR YTRANS=Y0 NASR ZTRANS=Z0 NASR PRINT592 NASR PRINT5003 NASR PRINT30,GELLIP,A,B,FINV,E NASR PRINT592 NASR PRINT40,XTRANS,YTRANS,ZTRANS NASR IF(CGELIP.EQ.ELNAME) GO TO 5074 NASR C NASR C PRINT THE PARAMETERS OF THE ARBITRARILY CHOZEN GEOCENTRIC ELLIPSOID NASR C NASR PRINT592 NASR PRINT5004 NASR PRINT30,CGELIP,ACG,BCG,FINVCG,ECG NASR C C 5074 CONTINUE NASR C NASR C PRINT THE APPROXIMATE GEOGRAPHICAL COORDINATES (LAT,LONG,HT) OF THE NASR C TRACKING STATION AND THE NECESSARY DESCRIPTION NASR C NASR PRINT592 NASR PRINT41 NASR PRINT42 NASR C NASR C RESET (ISTN = ICHOOZ) ,TO PRINT THE DATA FOR THE CHOZEN STATION ONLY NASR C NASR ISTN=ICHOOZ NASR PRINT43,STNAME,ISTN,IDEG1,MINT,SEC1,IDEG2,MIN2,SEC2,OBSHT NASR PRINT44,GELLIP NASR C NASR C PRINT THE APPROXIMATE CARTESIAN COORDINATES (X,Y,Z) OF THE TRACKING NASR C STATION WITH THEIR DESCRIPTION. NASR C NASR PRINT592 NASR PRINT45 NASR PRINT46 NASR PRINT47,STNAME,ISTN,OBSX,OBSY,OBSZ NASR PRINT49,GELLIP NASR C NASR C***********************************************************************NASR C NASR PRINT2000,KPASS,STNAME,ISTN NASR 2000 FORMAT('1'//,5X,'TOTAL NUMBER OF THE ACCEPTED PASSES ARE =',I5,3X,NASR 1'PASSES',3X,'AT *** ',A4,' *** STATION-WHOSE NUMBER IS :',I2/) NASR PRINT2073,IYEAR NASR 2073 FORMAT(40X,'YEAR OF OBSERVATIONS WAS : ',I4/) NASR C NASR C PRINT THE IDENTIFICATION CARDS AND THE MAIN CHARACTERISTICS (W.R.T. NASR C THE TRACKING STATION) OF ALL THE ACCEPTED PASSES ; ARRANGED INTO NASR C GROUPS OF FOUR (ONE FROM EACH QUADRANT) NASR C NASR PRINT2006 NASR 2006 FORMAT(4X,'NOTE : THE IDENTIFICATION CARD OF EACH PASS IN THE FOLLNASR 1OWING TABLE READS FROM RIGHT TO LEFT AS FOLLOWS : -'/) NASR PRINT2007 NASR 2007 FORMAT(10X,'3-DIGITS : DAY NUMBER'/10X,'2-DIGITS : LOCKON HOUR'/10NASR 1X,'2-DIGITS : LOCKON MINUTE'/10X,'2-DIGITS : SATELLITE NUMBER'/10XNASR 2,'2-DIGITS : PASS NUMBER W.R.T. THE DAY #'/10X,'1-DIGIT : NUMBER NASR 3OF ACCEPTED DOPPLERS'//) NASR PRINT2001 NASR 2001 FORMAT(30X,'CHARACTERISTICS OF ACCEPTED PASSES'/30X,'=============NASR 1====================='/) NASR PRINT2002 NASR 2002 FORMAT(10X,'SEQUENCE',5X,'IDENTIFICATION',6X,'MAXIMUM',8X,'TRACK',NASR 17X,'AZIMUTH',5X,'QUADRANT'/11X,'NUMBER',11X,'CARD',10X,'ELEVATION'NASR 2,5X,'DIRECTION'/) NASR JCOUNT=0 NASR DO 2003 JK=1,KPASS NASR PRINT2004,JK,IDL(JK),ISPE(JK),VAMSUM(JK),IDRSUM(JK),IAZSUM(JK),IQUNASR 1ADR(JK) NASR 2004 FORMAT(12X,I3,9X,I7,I5,8X,F5.2,10X,A3,10X,A4,8X,A3) NASR C C JCOUNT=JCOUNT+1 NASR IF(JCOUNT.NE.4) GO TO 2003 NASR PRINT2052 NASR JCOUNT=0 NASR 2003 CONTINUE NASR C NASR C***********************************************************************NASR C NASR C PRINT THE (XYZ) COORDINATES OF THE SATELLITE POSITIONS ,AND THE NASR C CORRESPONDING DOPPLERS ; FOR THE ACCEPTED PASSES. NASR C NASR PRINT2010,NAROWS,NACOLS NASR 2010 FORMAT('1'//10X,'NAROWS =',I5,2X,': IS THE TOTAL NUMBER OF ACCEPTENASR 1D DOPPLER EQUATIONS'//10X,'NACOLS =',I5,2X,': IS THE TOTAL NUMBER NASR 2OF UNKNOWN PARAMETERS'//) NASR PRINT2011 NASR 2011 FORMAT(5X,'TABLE OF SATELLITE COORDINATES (IN THE AVERAGE TERRESTRNASR 1IAL SYSTEM) AND THE CORRESPONDING DOPPLERS'/5X,'==================NASR 2==================================================================NASR 3============='/) NASR PRINT2012 NASR 2012 FORMAT(3X,'SEQUENCE NUMBER OF',3X,'SEQUENCE NUMBER OF',4X,'X -(METNASR 1ERS)',5X,'Y -(METERS)',5X,'Z -(METERS)',4X,'ACCEPTED DOPPLERS'/5X,NASR 2'SAT POSITIONS',6X,'ACCEPTED DOPPLERS') NASR IJ=1 NASR JJJ=0 NASR DO 2020 JJ=1,KPASS NASR PRINT2021,JJ NASR 2021 FORMAT(//30X,'*** P A S S S E Q U E N C E N U M B E R',I5,' ***'NASR 1/34X,'------------------------------------------'/) NASR IS=0 NASR JJJ=JJJ+NUMXYZ(JJ) NASR JJK=JJJ-1 NASR DO 2022 II=IJ,JJK NASR IS=IS+1 NASR PRINT2023,IS,(XSAT(II,J),J=1,3),IEQN(II+1),DOP(II+1) NASR 2023 FORMAT(10X,I3,28X,3(4X,F12.3)/31X,I1,63X,F11.2) NASR 2022 CONTINUE NASR ISJ=IS+1 NASR PRINT2024,ISJ,(XSAT(JJJ,J),J=1,3) NASR 2024 FORMAT(10X,I3,28X,3(4X,F12.3)) NASR IJ=JJJ+1 NASR 2020 CONTINUE NASR IREPET=IREPET+1 NASR IF(IREPET.LT.NCOPYS) GO TO 2030 NASR C NASR C***********************************************************************NASR C NASR C SCALING OF THE FREQUENCY OFFSET TO BE IN CYCLES PER 10 MINUTES ,FOR NASR C THE PURPOSE OF OBTAINING MORE CONSISTENT COEFFICIENTS OF THE DESIGN NASR C MATRIX-'A' ; NEEDED IN THE LEAST SQUARES ADJUSTMENT NASR C NASR DF010M=0.1920D+08 NASR C NASR C COMPUTE THE ACCUMULATED NUMBER OF OBSERVED SATELLITE POSITIONS FROM NASR C ALL ACCEPTED PASSES BEING PROCESSED HERE ; (NSUMJK) NASR C C NSUMJK=0 NASR DO 2090 ISUM=1,KPASS NASR NSUMJK=NSUMJK+NUMXYZ(ISUM) NASR 2090 CONTINUE NASR C NASR C***********************************************************************NASR C NASR C PUNCH AND PRINT THE NECESSARY DATA FOR THE (MATLAN)-PROGRAM 'BATCH' NASR C IF REQUIRED VIA THE OPTIONAL INDICATOR (IPUNCH) AS MENTIONED EARLIER NASR C NASR IF(IPUNCH.NE.1) GO TO 7777 NASR PRINT2050 NASR 2050 FORMAT('1'//27X,'PUNCHED DATA FOR THE (MATLAN)-PROGRAM ** B A T C NASR 1H **'/27X,'====================================================='/NASR 2) NASR C NASR C PUNCH AND PRINT THE TOTAL NUMBER OF ACCEPTED PASSES ,AND THE NASR C ACCUMULATED NUMBER OF OBSERVED SATELLITE POSITIONS FROM THESE PASSES NASR C NASR PUNCH2091,KPASS,NSUMJK NASR PRINT2091,KPASS,NSUMJK NASR 2091 FORMAT(2X,'NPASS =',I4/2X,'NSUMJK =',I4) NASR C NASR C PUNCH AND PRINT THE ROW AND COLUMN DIMENSION OF THE DESIGN MATRIX-A NASR C NASR PUNCH2051,NAROWS,NACOLS NASR PRINT2051,NAROWS,NACOLS NASR 2051 FORMAT(2X,'NAROWS =',I4/2X,'NACOLS =',I4) NASR C NASR C PUNCH AND PRINT THE YEAR OF OBSERVATIONS ,AND THE APRIORI VARIANCE- NASR C FACTOR (WHICH WAS ASSUMED IN CYCLES/SEC SQUARED UNITS) NASR C NASR PUNCH2072,IYEAR,IAVF NASR PRINT2072,IYEAR,IAVF NASR 2072 FORMAT(2X,'IYEAR =',I5/2X,'AVF =',I5) NASR PRINT2052 NASR 2052 FORMAT(//) NASR C NASR C***********************************************************************NASR C NASR C PUNCH AND PRINT THE WEIGHT COEFFICIENT MATRIX 'ELB' OF THE OBSERVED NASR C DOPPLER COUNTS ; AS AN IDENTITY MATRIX ( IDCORR = 0 ) NASR C NASR IF(IDCORR.EQ.1) GO TO 2070 NASR PUNCH2054,NAROWS,NAROWS,NAROWS,NAROWS NASR PRINT2054,NAROWS,NAROWS,NAROWS,NAROWS NASR 2054 FORMAT(2X,'ELB',I5,',',I3,2X,'N =',I4,2X,'D'/1X,'D*1,1*',2X,I3,'*0NASR 1.1D 01'/'END') NASR PRINT2052 NASR GO TO 5033 NASR 2070 CONTINUE NASR C NASR C RE-ARRANGE THE VECTOR (IEQN) WHOSE ELEMENTS INDICATE THE SEQUENCE NASR C NUMBER OF ACCEPTED DOPPLERS (LOCAL PER EACH PASS)WITH ZERO ELEMENTS NASR C FOR THE REJECTED DOPPLERS ; TO ANOTHER VECTOR (IEQNON) WHOSE ELEMENTSNASR C WILL BE THE NON-ZERO ELEMENTS OF IEQN-VECTOR. THIS IS DONE FOR THE NASR C C C PURPOSE OF INTRODUCING THE CORRELATION BETWEEN THE CONSECUTIVE NASR C DOPPLERS IN THE ELB-MATRIX NASR C NASR ICOUNT=0 NASR DO 2092 I=1,NSUMJK NASR IF(IEQN(I).EQ.0) GO TO 2092 NASR ICOUNT=ICOUNT+1 NASR IEQNON(ICOUNT)=IEQN(I) NASR 2092 CONTINUE NASR C NASR C ESTABLISH THE COFACTORS DIAGONAL ,WHICH WILL CONSTITUTE THE UPPER NASR C AND LOWER PARALLELS TO THE MAIN DIAGONAL OF ELB-MATRIX ,WITH ITS NASR C DIMENSION BEING LESS BY ONE THAN THE MAIN DIAGONAL SIZE (ICOUNT) NASR C -THIS IS DONE BY ASSIGNING A COFACTOR OF -0.5 FOR THE CONSECUTIVE NASR C DOPPLERS ,AND A COFACTOR OF 0 FOR THE NON-CONSECUTIVE ONES. NASR C NASR DO 2093 J=2,ICOUNT NASR COFACT(K)=0.0 NASR GO TO 2093 NASR 2099 COFACT(K)=-0.5 NASR 2093 CONTINUE NASR C NASR C PUNCH ,AND PRINT THE WEIGHT-COEFFICIENT MATRIX OF THE ACCEPTED NASR C DOPPLER COUNTS WITH ITS COFACTORS (OFF-DIAGONAL ELEMENTS) IN CASE NASR C OF CORRELATED ADJACENT DOPPLERS : IDCORR = 1 NASR C NASR C COMPUTE THE TOTAL NUMBER OF THE ELEMENTS FROM ELB-MATRIX ,NEEDED TO NASR C BE PUNCHED FOR THE MATLAN-PROGRAM 'BATCH' NASR C NASR NCOF=NAROWS-1 NASR NELB=NAROWS+2*NCOF NASR PUNCH2094,NAROWS,NAROWS,NELB,NAROWS NASR PRINT2094,NAROWS,NAROWS,NELB,NAROWS NASR 2094 FORMAT(2X,'ELB',I5,',',I3,2X,'N =',I4,2X,'D'/1X,'D*1,1*',2X,I3,'*0NASR 1.1D 01'/1X,'D*1,2*') NASR PUNCH2095,(COFACT(I),I=1,NCOF) NASR PRINT2095,(COFACT(I),I=1,NCOF) NASR 2095 FORMAT(8D10.1) NASR PUNCH2096 NASR PRINT2096 NASR 2096 FORMAT(1X,'D*2,1*') NASR PUNCH2095,(COFACT(I),I=1,NCOF) NASR PRINT2095,(COFACT(I),I=1,NCOF) NASR PUNCH2097 NASR PRINT2097 NASR 2097 FORMAT('END') NASR PRINT2052 NASR 5033 CONTINUE NASR C NASR C PUNCH AND PRINT THE APPROXIMATE VALUES OF THE UNKNOWN PARAMETERS ; NASR C (XYZ) OF THE TRACKING STATION AND THE FREQUENCY OFFSETS NASR C NASR PUNCH2053,NACOLS,NACOLS,KPASS,DF010M,OBSX,OBSY,OBSZ NAsR PRINT2053,NACOLS,NACOLS,KPASS,DF010M,OBSX,OBSY,OBSZ NAsR C C 2053 FORMAT(2X,'X0',I5,',1',2X,'N =',I4,2X,'D'/1X,'R*1,1*',2X,I3,'*',D1NASR 16.9,3(1X,D16.9)/'END') NASR C NASR C***********************************************************************NASR C NASR C PUNCH AND PRINT THE NUMBER OF OBSERVED SATELLITE POSITIONS PER EACH NASR C ACCEPTED PASS ,WITH THE RESULTING SATELLITE COORDINATES IN THE NASR C AVERAGE TERRESTRIAL SYSTEM ,AND THE CORRESPONDING DOPPLER COUNTS. NASR C NASR KK=0 NASR MM=1 NASR DO 2055 LL=1,KPASS NASR PUNCH2056,NUMXYZ(LL) NASR PRINT2056,NUMXYZ(LL) NASR 2056 FORMAT(2X,I3) NASR KK=KK+NUMXYZ(LL) NASR DO 2057 NN=MM,KK NASR PUNCH2058,(XSAT(NN,J),J=1,3),DOP(NN) NASR PRINT2058,(XSAT(NN,J),J=1,3),DOP(NN) NASR 2058 FORMAT(1X,D18.11,3(2X,D18.11)) NASR 2057 CONTINUE NASR MM=KK+1 NASR 2055 CONTINUE NASR C NASR C PUNCH AND PRINT THE MAIN CHARACTERISTICS OF THE ACCEPTED PASSES NASR C NASR PRINT2052 NASR DO 2060 KLM=1,KPASS NASR PUNCH2061,IDL(KLM),ISPE(KLM),VAMSUM(KLM),IDRSUM(KLM),IAZSUM(KLM),INASR 1QUADR(KLM) NASR PRINT2061,IDL(KLM),ISPE(KLM),VAMSUM(KLM),IDRSUM(KLM),IAZSUM(KLM),INASR 1QUADR(KLM) NASR 2061 FORMAT(10X,I7,I5,5X,F5.2,5X,A3,5X,A4,5X,A3) NASR 2060 CONTINUE NASR C NASR C PUNCH AND PRINT THE NAMES OF THE GEODETIC REFERENCE ELLIPSOID ,AND NASR C THE ARBITRARILY CHOZEN GEOCENTRIC ELLIPSOID ; RESPECTIVELY NASR C NASR PRINT2052 NASR PUNCH5079,GELLIP,CGELIP NASR PRINT5079,GELLIP,CGELIP NASR 5079 FORMAT(2(2X,A3)) NASR C NASR C PUNCH AND PRINT THE PARAMETERS OF THE GEODETIC ELLIPSOID ; SEMI- NASR C MAJOR AXIS ,RECIPROCAL OF FLATTENNING ,AND TRANSLATIONS FROM THE NASR C GEOCENTRE (WITH THE INDICATOR OF THEIR EXISTANCE) NASR C NASR PRINT2052 NASR PUNCH2065,A,FINV,XTRANS,YTRANS,ZTRANS,IDTRAN NASR PRINT2065,A,FINV,XTRANS,YTRANS,ZTRANS,IDTRAN NASR 2065 FORMAT(1X,F10.2,2X,F10.6,3(2X,F7.2),4X,I1) NASR IF(CGELIP.EQ.ELNAME) GO TO 5077 NASR C C PUNCH AND PRINT THE SEMI-MAJOR AXIS AND THE RECIPROCAL OF FLATTENNINGNASR C OF THE CHOZEN GEOCENTRIC ELLIPSOID (IF ANY) ; NASR C NASR C C C PRINT2052 NASR PUNCH5078,ACG,FINVCG NASR PRINT5078,ACG,FINVCG NASR 5078 FORMAT(1X,F10.2,2X,F10.6) NASR 5077 CONTINUE NASR C NASR C PUNCH AND PRINT THE NAME & NUMBER OF THE CHOZEN TRACKING STATION NASR C NASR PRINT2052 NASR PUNCH2066,STNAME,ISTN NASR PRINT2066,STNAME,ISTN NASR 2066 FORMAT(1X,A4,2X,I2) NASR 7777 CONTINUE NASR STOP NASR END NASR C C C C*************************************************************************** C* * C* FUNCTION 'RADIAN' : * C* ----------------- * C* THE PURPOSE OF THIS FUNCTION SUBROUTINE IS TO * C* CONVERT THE GIVEN ANGLE IN (DEG,MIN & SEC) INTO ITS EQUIVALENT VALUE * C* IN RADIANS TO BE USED IN OTHER COMPUTATION. * C* * C* THIS FUNCTION SUBROUTINE HAS AN ENTRY-POINT 'IDEG' TO PERFORM THE * C* INVERSE OF THE ABOVE CASE , NAMELY TO CONVERT THE GIVEN ANGLE FROM * C* RADIANS TO ITS EQUIVALENT VALUE IN (DEG,MIN & SEC). * C* * C* NOTICE : IF THE INPUTTED ANGLE IN RADIANS IS NEGATIVE ,THE -VE SIGN * C* WILL APPEAR ONLY BESIDE THE DEGREES AFTER CONVERSION ,AND THIS WILL * C* IMPLY THAT THE WHOLE ANGLE (DEG,MIN & SEC) IS -VE. * C* * C* THIS FUNCTION SUBROUTINE AND ITS ENTRY-POINT ARE VALID FOR ONLY ONE * C* CONVERSION AT A TIME. * C* * C*************************************************************************** C FUNCTION RADIAN(DEG,XMIN,SEC) IMPLICIT REAL*8(A-H,O-Z) PI=3.141592653589793 RADIAN=(DEG+XMIN/60.+SEC/3600.)*PI/180. RETURN C C*************************************************************************** C ENTRY IDEG(RAD,MIN,SEC) PI=3.141592653589793 SIGN=RAD RAD=DABS(RAD) DEG=RAD*180./PI+0.005/3600. IDEG=DEG MIN=(DEG-IDEG)*60. SEC=((DEG-IDEG)*60.-MIN)*60.-0.005 IF(SIGN.LT.0.) IDEG=-IDEG AMIN=MIN AMIN=DABS(AMIN) MIN=AMIN SEC=DABS(SEC) RAD=SIGN RETURN END C C C C C*************************************************************************** C* * C* SUBROUTINE 'GEOREC' : * C* ------------------- * C* THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C* (X,Y,Z) GEOCENTRIC CARTESIAN COORDINATES OF A SINGLE STATION ,REFERRED * C* TO THE CENTRE OF GRAVITY OF THE EARTH AS THE ORIGIN OF THE COORDINATE * C* SYSTEM ; FROM ITS GIVEN GEODETIC COORDINATES (LAT,LONG & HT) WHICH ARE * C* REFERRING TO A GIVEN REFERENCE ELLIPSOID. * C* * C* THE REQUIRED INPUT DATA ARE : * C* ---------------------------- * C* * C* 1-THE PARAMETERS OF THE REFERENCE ELLIPSOID ; SEMI-MAJOR AXIS (A) & * C* SEMI-MINOR AXIS (B) ,IN METERS. * C* 2-THE TRANSLATION COMPONENTS (X0,Y0,Z0) OF THE REFERENCE ELLIPSOID * C* FROM THE GEOCENTRE -(IN METERS) * C* 3-THE GEODETIC (LAT & LONG) IN RADIANS * C* 4-THE GEOMETRICAL HEIGHT ABOVE THE REFERENCE ELLIPSOID -(IN METERS) * C* * C* THE OUTPUT WILL BE THE (X,Y,Z) GEOCENTRIC CARTESIAN COORDINATES OF THE * C* GIVEN STATION - IN METERS. * C* * C* NOTE : THIS SUBROUTINE IS VALID FOR ONLY ONE STATION AT A TIME * C* * C*************************************************************************** C SUBROUTINE GEOREC(A,B,XO,YO,ZO,PHI,XLAMB,HT,X,Y,Z) IMPLICIT REAL*8(A-H,O-Z) BOA=(B/A)*(B/A) XN=A/DSQRT(DCOS(PHI)*DCOS(PHI)+BOA*DSIN(PHI)*DSIN(PHI)) X=XO+(XN+HT)*DCOS(PHI)*DCOS(XLAMB) Y=YO+(XN+HT)*DCOS(PHI)*DSIN(XLAMB) Z=ZO+(XN*BOA+HT)*DSIN(PHI) RETURN END C C C*************************************************************************** C* * C* SUBROUTINE 'RECGEO' : * C* ------------------- * C* THE PURPOSE OF THIS SUBROUTINE IS TO CONVERT THE * C* (X,Y,Z) GEOCENTRIC CARTESIAN COORDINATES OF A SINGLE STATION ,TO ITS * C* CORRESPONDING (LAT,LONG,HT) W.R.T. A GIVEN REFERENCE ELLIPSOID (EITHER* C* GEODETIC OR GEOCENTRIC) ; USING THE ITERATIVE PROCEDURE. * C* * C* THE REQUIRED INPUT DATA ARE : * C* ---------------------------- * C* 1-THE (X,Y,Z) GEOCENTRIC CARTESIAN COORDINATES - (IN METERS) * C* 2-THE SEMI-MAJOR & SEMI-MINOR AXES (A & B) OF THE REFERENCE ELLIPSOID * C* -(IN METERS) * C* 3-THE TRANSLATION COMPONENTS (X0,Y0,Z0) OF THE REFERENCE ELLIPSOID * C* FROM THE GEOCENTRE -(IN METERS) * C* * C* THE OUTPUT WILL BE THE LATITUDE (PHI) & LONGITUDE (XLAMB) -IN RADIANS * C* AND THE GEOMETRICAL HEIGHT (H) ABOVE THE ELLIPSOID -IN METERS ,REFERR- * C* ING TO THE GIVEN REFERENCE ELLIPSOID. * C* * C* NOTE : THIS SUBROUTINE IS VALID ONLY FOR ONE STATION AT A TIME. * C*************************************************************************** C SUBROUTINE RECGEO(A,B,XO,YO,ZO,X,Y,Z,PHI,XLAMB,HT) IMPLICIT REAL*8(A-H,O-Z) X1=X-XO Y1=Y-YO Z1=Z-ZO BOA=(B/A)*(B/A) E2=1.-BOA P=DSQRT(X1*X1+Y1*Y1) XLAMB=DATAN2(Y1,X1) C C INITIAL APPROXIMATIONS FOR N & H C XN=A HT=DSQRT(X1*X1+Y1*Y1+Z1*Z1)-A PHI=DATAN((Z1/P)/(1.-E2*XN/(XN+HT))) C C ITERATIVE LOOP TO CONVERGE TO ONE MILLIMETER. C ICOUNT=0 100 CONTINUE ICOUNT=ICOUNT+1 HT1=HT PHI1=PHI XN=A/DSQRT(DCOS(PHI)*DCOS(PHI)+BOA*DSIN(PHI)*DSIN(PHI)) HT=P/DCOS(PHI)-XN PHI=DATAN((Z1/P)/(1.-E2*XN/(XN+HT))) IF(ICOUNT.GT.20) GOTO 200 IF(DABS(HT1-HT).GT.A/1E10) GOTO 100 IF(DABS(PHI1-PHI).GT.1E-10) GOTO 100 200 CONTINUE RETURN END C C C C*************************************************************************** C* * C* SUBROUTINE 'DROTAT' : * C* ------------------- * C* THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE * C* (3,3) PRODUCT ROTATION MATRIX (ROT) ,RESULTING FROM A SERIES OF UP TO * C* NINE CONSECUTIVE ROTATIONS ,REQUIRED FOR THE COORDINATE TRANSFORMATION * C* BETWEEN THE DIFFERENT SYSTEMS. * C* * C* THE REQUIRED INPUT DATA ARE : * C* ---------------------------- * C* * C* 1-ANGLE : IS A VECTOR OF UP TO 9 CONSECUTIVE ANGLES OF SINGLE AXIS * C* ROTATIONS - (IN RADIANS) * C* 2-NAXIS : IS A VECTOR OF THE AXIS OF ROTATION NUMBERS (EITHER 1 , 2 OR * C* 3 ,FOR EACH GIVEN ANGLE) ,I.E. SAME SIZE AS THE VECTOR-ANGLE * C* 3-NUM : IS THE NUMBER OF INPUTTED ROTATIONS ,WHICH IS THE SIZE OF * C* BOTH (ANGLE & NAXIS)-VECTORS. * C* * C* THE OUTPUT WILL BE THE (3,3) PRODUCT ROTATION MATRIX (ROT) ,COMPUTED * C* IN DOUBLE PRECISION * C* * C* NOTE : THIS SUBROUTINE IS VALID UP TO 9 CONSECUTIVE SINGLE ROTATIONS * C* * C*************************************************************************** C SUBROUTINE DROTAT(ANGLE,NAXIS,NUM,ROT) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ROT(3,3),R1(3,3),R2(3,3),ANGLE(NUM),NAXIS(NUM) DO 100 K=1,NUM C C FIND OTHER TWO AXES N2,N3 C N2=NAXIS(K)+1 IF(N2.GT.3) N2=N2-3 N3=N2+1 IF(N3.GT.3) N3=N3-3 C C SET DIAGONAL ELEMENTS C R1(NAXIS(K),NAXIS(K))=1. R1(N2,N2)=DCOS(ANGLE(K)) R1(N3,N3)=R1(N2,N2) C C SET NONZERO OFF-DIAGONAL ELEMENTS C R1(N2,N3)=DSIN(ANGLE(K)) R1(N3,N2)=-R1(N2,N3) C C SET ZERO OFF-DIAGONAL ELEMENTS C R1(NAXIS(K),N2)=0. R1(NAXIS(K),N3)=0. R1(N2,NAXIS(K))=0. R1(N3,NAXIS(K))=0. C C C C IF FIRST ROTATION, SET ROT=R1 C IF(K.GT.1) GOTO 20 DO 10 I=1,3 DO 10 J=1,3 10 ROT(I,J)=R1(I,J) GOTO 100 C C IF NOT FIRST ROTATION SET R2=R1*ROT AND ROT=R2 C 20 DO 30 J=1,3 DO 30 I=1,3 R2(I,J)=0. DO 30 M=1,3 30 R2(I,J)=R2(I,J)+R1(I,M)*ROT(M,J) DO 40 I=1,3 DO 40 J=1,3 40 ROT(I,J)=R2(I,J) 100 CONTINUE RETURN END C C C C*************************************************************************** C* * C* FUNCTION FOR READING FROM CARDS FORMAT FREE * C* * C* R.L.PARKER (SEE BEDFORD INSTITUTE COMPUTER * C* NOTE 1969-8-C FOR DETAILS) * C* * C* MODIFICATIONS BY D.E.WELLS UNB JAN 1971 * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C**********PARAMETER DESCRIPTIONS********** * C* * C* INPUT IN =0 NO PRINT * C* =1 PRINT MESSAGE AND BAD NUMBERS ONLY * C* =2 PRINT ALL * C* =3 REREAD LAST NUMBER AND PRINT * C* * C* OUTPUT IOUT=0 IF SUCCESSFUL READ * C* = FIRST BAD CHARACTER IF ERROR HIT * C* * C* COLUMN COUNTER ICOL AVAILABLE AS ARGUMENT * C* USUALLY LEFT UNTOUCHED BY USER * C* HOWEVER FOR FLEXIBILITY IT CAN BE SET WHEN FUNCTION CAL * C* ICOL=81 WILL SKIP TO END OF CURRENT CARD * C* ICOL.LT.81 WILL SKIP BACK OR FORWARD TO THAT COLUMN * C* * C**********POSSIBLE MODIFICATIONS********** * C* * C* IF 'IN' NOT REQUIRED - * C* DROP FROM PARAMETER LISTS IN FUNCTION STATEMENT AND IN ALL * C* CALLING STATEMENTS * C* INITIALIZE TO DESIRED VALUE ( 0 1 2 OR 3) USING EITHER * C* DATA STATEMENT ( E.G. 'DATA IN/1/ ' ) OR * C* REPLACEMENT STATEMENT (E.G. 'IN=2 ' ) * C* * C* IF ' IOUT' NOT REQUIRED - * C* DROP FROM PARAMETER LISTS IN FUNCTION STATEMENT AND IN ALL * C* CALLING STATEMENTS * C* * C* IF 'ICOL' NOT REQUIRED - * C* DROP FROM PARAMETER LISTS IN FUNCTION STATEMENT AND IN ALL * C* CALLING STATEMENTS * C* INITIALIZE TO 81 USING DATA STATEMENT ( DATA ICOL/81/ ) * C* * C* IF DOUBLE PRECISION DATA EXPECTED (MORE THAN SEVEN DIGITS) - * C* ADD TO BOTH FUNCTION AND TO CALLING PROGRAM THE STATEMENT * C* ' DOUBLE PRECISION FREAD ' * C* * C*************************************************************************** C FUNCTION FREAD(IN,IOUT) DOUBLE PRECISION X,XX,B DOUBLE PRECISION FREAD DIMENSION IA(80),IDIGIT(17) C C DATA IDIGIT/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9, 1 1H ,1H, ,1H.,1HE,1HD,1H-,1H+/ DATA IE/3/ DATA ICOL/81/ C*****TEST FOR REREAD MODE (IN=3) IF(IN.GE.3) GOTO(400,410,860),IE IOUT=0 C*****TEST FOR EMPTY CARD BUFFER. READ NEW CARD IF EMPTY IF(ICOL.GT.80) GOTO 201 C** IE IS ONE IF NO E YET HIT C IC1 ADDRESSES FIRST CHARACTER IN FREAD CALL C B IS POWER OF TEN BY WHICH FRACTIONAL PART MUST BE DIVIDED 100 IE=1 IC1=ICOL B=1.0 C*****RETURN HERE TO PICK UP EXPONENT C IDEC IS ONE IF NO DECIMAL OR E YET HIT C X IS NUMBER EVALUATED (FREAD) C SIGN IS SIGN OF FREAD C NC IS LAST ELEMENT ON IDIGIT LIST 110 IDEC=1 X=0.0 SIGN=+1.0 NC=17 C** IGNORE LEADING SPACES AND COMMAS DO 150 IC=ICOL,80 IF(IA(IC).NE.IDIGIT(11).AND.IA(IC).NE.IDIGIT(12)) GOTO 202 150 CONTINUE C*****END OF CARD. ERROR IF EXPONENTIAL PART EXPECTED. ICOL=ICOL-1 IF(IE.GT.1) GOTO 800 C** NO ERROR. READ ANOTHER CARD 201 READ 2000,(IA(J),J=1,80) 2000 FORMAT(80A1) ICOL=1 GO TO 100 C **************************************************************** C NO MORE LEADING SPACES AND COMMAS C START TESTING CHARACTER BY CHARACTER UNTIL TERMINATING SPACE C OR COMMA HIT 202 DO 290 ICOL=IC,80 JA=IA(ICOL) C COMPARE EACH CHARACTER TO IDIGIT LIST DO 205 I=1,NC IF(JA.EQ.IDIGIT(I)) GOTO 206 205 CONTINUE C CHARACTER NOT ON LIST. ERROR GO TO 800 C CHARACTER ON LIST. TEST WHETHER NUMERAL 206 J=I-9 IF(J.LE.0) GOTO 25 GOTO(25,300,300,210,310,310,220,290),J C CHARACTER IS DECIMAL. ERROR IF PREVIOUS DECIMAL OR E HIT 210 IDEC=IDEC+IE IF(IDEC.GT.2) GOTO 800 GOTO 290 C C C CHARACTER IS MINUS. RESET SIGN AND SHORTEN IDIGIT LIST 220 SIGN=-1.0 GO TO 290 C CHARACTER IS NUMERAL 25 DIGIT=I-1 GOTO(251,252,800),IDEC C INTEGRAL PART OF NUMBER 251 X=10.0D0*X+DIGIT GO TO 290 C FRACTIONAL PART OF NUMBER 252 B=B/10.0D0 X=X+B*DIGIT C CHARACTER IS PLUS. SHORTEN IDIGIT LIST 290 NC=15 C ***************************************************************** C NO TERMINATING SPACE OR COMMA BEFORE END OF CARD HIT ICOL=81 C CHARACTER IS TERMINATING SPACE OR COMMA. TEST FOR E OR F FORMAT 300 GO TO(400,410),IE C CHARACTER IS E. SET XX=MANTISSA 310 IE=IE+1 XX=X*SIGN C ERROR IF NO MANTISSA IF(ICOL.EQ.IC) GOTO 800 ICOL=ICOL+3-IE GO TO (110,110,800),IE C TERMINATION OF F FORMAT 400 FREAD=X*SIGN GO TO 500 C TERMINATION OF E FORMAT 410 IX=X*SIGN FREAD=XX*10.0D0**IX C** PRINT IF IN=2 OR 3 500 IC2=ICOL-1 IF (IN .GT.1) PRINT 5000,(IA(J),J=IC1,IC2) 5000 FORMAT (10X,80A1) RETURN C ***************************************************************** C ERROR HANDLING 800 IOUT=IA(ICOL) C SEARCH FOR NEXT TERMINATING BLANK OR COMMA DO 810 I=ICOL,80 IF(IA(I).EQ.IDIGIT(11).OR.IA(I).EQ.IDIGIT(12)) GO TO 815 810 CONTINUE I=81 815 ICOL=I IF(IN.EQ.0) GOTO 850 IC2=ICOL-1 PRINT 820,IOUT,(IA(J),J=IC1,IC2) 820 FORMAT(25H MALFORMED NUMBER DUE TO ,A1,4H IN ,80A1) 850 GO TO(400,410,400),IE 860 PRINT 870 870 FORMAT(25H ILLEGAL REREAD ATTEMPTED ) FREAD=0. RETURN END