C PROGRAM READSEDR C 05/02/97: user notes C This FORTRAN code was originally named SEDREXT and was supplied C by Ruedi von Steiger of the SWICS team at the University of Bern. C See the SEDR file documentation for definitions of word contents. C The code has been validated at NSSDC for VAX operation under VMS. C C Input file: SEDRyymm.FIX {enter only the prefix when prompted} for C year yy and month mm. These files consist of binary data C records with data fields in binary format. C Start record: enter "1" to start at the first input record. C Record frequency: enter "1" to get all records. C SEDR word number list: enter up to 32 word numbers from the SEDR C document Table 4. Numbers 163-168 give the six C elements of the state vector (3 position + 3 C velocity components) for the Ulysses spacecraft C in Jupiter-centered System III coordinates. C Output file: SEDRyymm.INF (year and fractional day of year + ASCII list C format for selected SEDR words. The FORTRAN output format C should be changed to suit the needs of the user. C Original documentation: C program to extract parameters as a C function of sc event time form ulysses sedr tapes. C written by r.v.steiger, adapted from fsedr5.for C july 22, 1992 C-------------------------- DECLARATION PART --------------------------- IMPLICIT NONE C each supplemetary experiment data record consists of 1152 bytes C grouped into 288 integer*4 variables INTEGER*4 LSEDR PARAMETER( LSEDR=1152 ) INTEGER*2 I2SEDR(LSEDR/2) INTEGER*4 I4SEDR(LSEDR/4) EQUIVALENCE (I2SEDR,I4SEDR) CHARACTER*1152 CSEDR EQUIVALENCE (I2SEDR,CSEDR) C buffer for reading from sedr file CHARACTER*512 BUFFER, TEMP C flag set while sedr is assembled from file LOGICAL*1 SEDRFULL /.FALSE./ C position where sedr ends in buffer, index INTEGER*2 SEDREND, SEDRI INTEGER*4 WORD(0:32) C first sedr possibly contains factors instead of data if facflg set INTEGER*4 FACTOR(26:288) LOGICAL*1 FACFLG C sedr parameters in real format, conversion routine to reals REAL RSEDR(26:288), IBMREAL C conversion factor from km to au REAL KMPAU PARAMETER( KMPAU=1.495979E8 ) C integer and fractional part of spacecraft event time INTEGER*4 ISCET, FSCET C floating point form of scet REAL*8 RSCET C function to convert time to floating point format REAL*8 RTIME EXTERNAL RTIME C time conversion routine and variables EXTERNAL SCET2DT INTEGER YY, MO, DD, HH, MI, SS, MS, DOY C directory and name of file CHARACTER*40 FILNAM C counters for valid and missing sedrs, output freq., dummy counter INTEGER*4 SEDRNR, MISSNR, FIRSTREC, NTHREC, I C routine to reverse byte order in integer*4 variables EXTERNAL REVBYT C--------------------------- STATEMENT PART ---------------------------- C print version info WRITE(*,*) ' SEDR decoding program, Version 1.6' WRITE(*,*) ' by r.v.steiger' WRITE(*,*) ' July 22, 1992' C initialize counters and flags SEDRNR = 0 MISSNR = 0 FACFLG = .FALSE. WRITE(*,2) ' ENTER NAME OF SEDR FILE: ' 02 FORMAT( A, $ ) ACCEPT 03, FILNAM 03 FORMAT( A ) WRITE(*,2) ' ENTER NUMBER OF SEDR WORDS TO EXTRACT (UP TO 32): ' ACCEPT 04, WORD(0) WRITE(*,2) ' ENTER THEIR NUMBERS: ' READ(*,*) ( WORD(I), I=1, WORD(0) ) WRITE(*,'(A,$)') ' ENTER NUMBER OF FIRST RECORD TO OUTPUT: ' ACCEPT 04, FIRSTREC WRITE(*,'(A,$)') ' ENTER OUTPUT FREQUENCY OF RECORDS: ' ACCEPT 04, NTHREC 04 FORMAT( I ) WRITE(*,05) FILNAM, FILNAM 05 FORMAT( ' READING SUPPLEMENTARY DATA RECORDS FROM ',A,'.SEDR',/, 2 ' WRITING SELECTED SEDR DATA TO ', A, '.INF' ) WRITE(*,*) C open files OPEN ( 10, FILE=' '//FILNAM//'.FIX', STATUS='OLD', READONLY ) OPEN ( 20, FILE=' '//FILNAM//'.INF', STATUS='NEW', 1 RECL=MAX(80,16+12*WORD(0)) ) C write header to output file WRITE(20,07) FILNAM, NTHREC 07 FORMAT( 'SUPPLEMENTARY DATA RECORDS FROM ', A, '.SEDR', 1 /, 'ONE RECORD PRINTED OUT OF', I4, / ) WRITE(20,08) 'YEAR', 'DOY', (WORD(I), I=1, WORD(0) ) 08 FORMAT( A4, X, A7, 32I12 ) C loop over all sedr on the tape DO WHILE (.TRUE.) C read sedr from tape (1152 bytes), exit if no more records available SEDRI = SEDREND IF ( SEDREND.NE.0 ) THEN CSEDR(1:SEDREND) = TEMP(1:SEDREND) END IF DO WHILE ( .NOT.SEDRFULL ) READ(10,10,END=9999) BUFFER 10 FORMAT( A512 ) DO I = 1, 512 SEDRI = SEDRI + 1 IF ( SEDRI.LE.LSEDR ) THEN CSEDR(SEDRI:SEDRI) = BUFFER(I:I) ELSE IF ( .NOT.SEDRFULL ) SEDREND = 0 SEDRFULL = .TRUE. SEDREND = SEDREND + 1 TEMP(SEDREND:SEDREND) = BUFFER(I:I) END IF END DO END DO SEDRFULL = .FALSE. C reorder bytes because tape was written on a ibm CALL REVBYT (I4SEDR, LSEDR/4 ) C extract spacecraft event time from header CALL MVBITS( I4SEDR(3), 0, 16, ISCET, 16 ) CALL MVBITS( I4SEDR(4), 16, 16, ISCET, 0 ) CALL MVBITS( I4SEDR(4), 0, 16, FSCET, 0 ) RSCET = RTIME( ISCET, FSCET ) CALL SCET2DT( RSCET, YY, MO, DD, HH, MI, SS, MS ) C calculate parameters from float format and store them in file IF (I4SEDR(24).EQ.0) THEN SEDRNR = SEDRNR + 1 WRITE(*,11) ' READING SEDR # ', SEDRNR, ' AT SCET', RSCET 11 FORMAT( A, I5, A, F15.3 ) DO I = 26, 288 RSEDR(I) = IBMREAL( I4SEDR(I) ) END DO IF ( SEDRNR.GE.FIRSTREC .AND. . MOD( SEDRNR-FIRSTREC, NTHREC ).EQ.0 ) THEN C write parameters to general output file WRITE(20,15) YY, DOY( YY, MO, DD ) + . HH/24E0+MI/1.440E3+SS/8.64E5+MS/8.64E8, . (RSEDR(WORD(I)), I=1, WORD(0)) 15 FORMAT( I4, 1X, F7.3, 1P, 32 E12.4 ) END IF C store scale factors or multiply by them, according to sedr-id ELSE IF (I4SEDR(24).EQ.1) THEN DO I = 26, 288 FACTOR(I) = I4SEDR(I) END DO FACFLG = .TRUE. C scale integers to proper value and write parameters to file ELSE IF ( I4SEDR(24).EQ.2 .AND. FACFLG ) THEN SEDRNR = SEDRNR + 1 IF ( MOD( SEDRNR-1, NTHREC ) .EQ. 0 ) THEN WRITE(*,20) ' READING SEDR # ', SEDRNR, ' AT SCET', RSCET 20 FORMAT( '+', A, I5, A, F15.3 ) END IF DO I = 26, 288 RSEDR(I) = I4SEDR(I) * 10E0**FACTOR(I) END DO IF ( SEDRNR.GE.FIRSTREC .AND. . MOD( SEDRNR-FIRSTREC, NTHREC ).EQ.0 ) THEN C write parameters to general output file WRITE(20,15) YY, DOY( YY, MO, DD ) + . HH/24E0+MI/1.440E3+SS/8.64E5+MS/8.64E8, . (RSEDR(WORD(I)), I=1, WORD(0)) END IF ELSE MISSNR = MISSNR + 1 END IF END DO 9999 CONTINUE C close files CLOSE( 10 ) CLOSE( 20 ) WRITE(*,40) SEDRNR, MISSNR 40 FORMAT( ' ', I5, ' SUPPLEMENTARY DATA RECORDS READ FROM TAPE', /, 1 ' ', I5, ' BLOCKS MISSING OR INVALID' ) STOP END C C C INTEGER*4 FUNCTION DOY( YY, MO, DD ) C C function to calculate the day number of a year given the date. C written by r.v.steiger C march 28, 1991 C IMPLICIT NONE C C input: year, month, day INTEGER*4 YY, MO, DD C C factor to distinguish normal and leap years INTEGER*4 ONETWO C C IF ( MOD(YY,4).EQ.0 .AND. MOD(YY,400).NE.0 ) THEN ONETWO = 1 ELSE ONETWO = 2 END IF C DOY = INT( 275*MO/9 ) - ONETWO * INT( (MO+9)/12 ) + DD - 30 C RETURN END C C C DOUBLE PRECISION FUNCTION RTIME( ITIME, FTIME ) C C function to convert the 6-byte time format to floating point C written by r.v.steiger C november 9, 1989; simplified october 28, 1991 C------------------------- DECLARATION PART ---------------------------- C integer and fractional part of time, dummy counter INTEGER*4 ITIME INTEGER*2 FTIME C--------------------------- STATEMENT PART ---------------------------- RTIME = 1D0 * ITIME + FTIME / 65536D0 RETURN END C C C SUBROUTINE SCET2DT( RSCET, YY, MO, DD, HH, MI, SS, MS ) C C subroutine to calculate the date and time from Ulysses SCET, i.e. C seconds since January 1, 1950 C written by r.v.steiger C march 26, 1991 C IMPLICIT NONE C C input: RSCET, Ulysses time in seconds since January 1, 1950 REAL*8 RSCET C C time in julian centuries since Jan. 0.5, 1900 REAL*8 T C C output: year, month, day, hour, minute, second, milliseconds INTEGER*4 YY, MO, DD, HH, MI, SS, MS C C T = RSCET / ( 36525D0 * 86400D0 ) + 5D-1 C CALL JUL2DT( T, YY, MO, DD, HH, MI, SS, MS ) C RETURN END C C C REAL FUNCTION IBMREAL( IN ) C C function to decode ibm-coded reals C written by r.v.steiger, oct. 28, 1991 C IMPLICIT NONE INTEGER*4 IN, MANTISSA, EXPONENT, SIGN REAL*4 R, MAXREAL /1.7014117E38/ C MANTISSA = IBITS( IN, 0, 24 ) EXPONENT = IBITS( IN, 24, 7 ) SIGN = IBITS( IN, 31, 1 ) C IF ( MANTISSA.LE.0 ) THEN R = 0E0 ELSE IF ( ( LOG(REAL(MANTISSA)) + (EXPONENT-70E0)*LOG(16E0) ) . .LT. LOG(MAXREAL) ) THEN R = (1-2*SIGN) * MANTISSA * 16E0**(EXPONENT-70) ELSE R = (1-2*SIGN) * MAXREAL TYPE *, ' OVERFLOW ERROR IN IBMREAL - RESULT SET TO MAXREAL' END IF C IBMREAL = R C RETURN END C C C SUBROUTINE REVBYT( IBUF, NUM ) C subroutine to reverse the order of bytes in an integer*4 array C ibuf of length num. C written by r.v.steiger C february 23, 1990 IMPLICIT NONE C input array, number of elements INTEGER*4 NUM, IBUF(NUM) C dummy variable, counters INTEGER*4 IDUM, I, J DO I = 1, NUM DO J = 0, 24, 8 CALL MVBITS( IBUF(I), J, 8, IDUM, 24-J ) END DO IBUF(I) = IDUM END DO RETURN END C C C SUBROUTINE JUL2DT( T, YY, MO, DD, HH, MI, SS, MS ) C C subroutine to calculate the date and time from julian centuries C since January 0.5, 1900. C For Ulysses, remember that T = SCET / (36525*86400) + 0.5 C written by r.v.steiger (using "astronomical formulae" by j.meeus) C march 5, 1990 C corrected march 26, 1991: specify all real constants in D format, C rather than E format. (rvst) C corrected again july 23, 1991: specify _really_ all real constants C in D format (rvst) C corrected november 1994: round instead of truncate the time to C integer values of milliseconds (R. Bodmer) C IMPLICIT NONE C C input: time in julian centuries since january 0.5, 1900; C julian day number REAL*8 T, JD C fractional and integer part of jd+0.5 INTEGER*4 Z REAL*8 F C C output: year, month, day, hour, minute, second, milliseconds INTEGER*4 YY, MO, DD, HH, MI, SS, MS C C auxiliary variables INTEGER*4 ALPHA, A, B, C, D, E C C calculate julian day number JD = T * 36525D0 + 2415020D0 JD = JD + 0.5D0 / 86400000D0 ! add half a millisecond for ! correct rounding (R. Bodmer, Nov. 94) Z = INT( JD + 5D-1 ) F = JD + 5D-1 - Z C C correct for gregorian calendar after oct. 4, 1582 IF ( Z .GE. 2299161 ) THEN ALPHA = INT( ( Z - 1867216.25D0 ) / 36524.25D0 ) A = Z + 1 + ALPHA - INT( ALPHA / 4D0 ) ELSE A = Z END IF C C calculate auxiliary results B = A + 1524 C = INT( ( B - 122.1D0 ) / 365.25D0 ) D = INT( 365.25D0 * C ) E = INT( ( B - D ) / 30.6001D0 ) C C find day, month, and year DD = B - D - INT( 30.6001D0 * E ) IF ( E .LT. 13.5 ) THEN MO = E - 1 ELSE MO = E - 13 END IF IF ( MO .GT. 2.5 ) THEN YY = C - 4716 ELSE YY = C - 4715 END IF C C find hour, minute, second, and millisecond HH = INT( 24 * F ) MI = INT( 1440 * F ) - 60 * HH SS = INT( 86400 * F ) - 60 * MI - 3600 * HH MS = INT( 86400000 * F ) - 1000 * SS - 60000 * MI - 3600000 * HH C RETURN END