	Options /Extend_Source
C**********************************************************************		
C       PROGRAM TO READ THE HAWKEYE MASTER SCIENCE FILE
C**********************************************************************
	common/trecord/id
	dimension bb0(9),bb1(9)
	character*80 filename,filename2,finter
        INTEGER*4 ID,IRAW
	integer*2 dummy 
	real*8 epoch0,epoch,ep,ep0,dep,tod,tstart,tave
C**********************************************************************
        REAL*8 TRJ,CD,EW,WORD24,WORD25,Z,ZZ,XI,YY,YI,Y1,Y2,Y3,EM,EI
C**********************************************************************
	dimension spde(16),spdb(8),freqs(16),dfreqs(16)
        DIMENSION ID(324),IRAW(64),TRJ(77),CD(21),EW(2),WORD24(16)
        DIMENSION WORD25(32),BX(3,32),BY(3,32),BZ(3,32),BCD(9,32)
        DIMENSION V(10),SA1(8,16),SA2(8,16),SFR(8,32),VD(6,64),IV(7)
        DIMENSION GM(32,2),FP(32,2),FE(32,2),ST(32,2),PP(8),RAMPE(16,16)
        DIMENSION RAMPP(16,16),IW24(16),IW25(32),EL(16,2),BB(3,32)
        DIMENSION SFF(8),SAF1(8),SAF2(8),AN(32,2),ANE(16,16),ANP(16,16)
        DIMENSION RPEN(16),RPEW(16),RPGF(16),ENLP(32),EWLP(32),ENPL(16)
        DIMENSION REEN(16),REEW(16),REGF(16),GFP(32),GFE(32),SFFBW(8)
        DIMENSION SA1BW(8),SA2BW(8),SA1C(8),SA2C(8),ANR(16),EM(9),EI(9)
C**********************************************************************
        EQUIVALENCE (ID(34),IRAW(1)),(ID(104),EW(1)),(ID(130),CD(1))
        EQUIVALENCE (ID(180),TRJ(1)),(ID(242),EM(1)),(ID(260),EI(1))
C************************************************************************
C       FREQUENCIES FOR THE VLF CHANNELS
        DATA SFF/13.3,17.8,23.7,31.1,42.2,56.2,100.,178./       !SFR KHZ
        DATA SAF1/1.78,5.62,17.8,56.2,1.78,5.62,17.8,56.2/      !SA-1 HZ
        DATA SAF2/178.,562.,1780.,5620.,178.,562.,1780.,5620./  !SA-2 HZ
C************************************************************************
C       BAND WIDTHS FOR VLF CHANNELS
        DATA SFFBW/878.,987.,1557.,1162.,1559.,2259.,4929.,5609./
        DATA SA1BW/2.80,1.51,13.25,26.46,3.13,1.98,13.28,24.08/
        DATA SA2BW/69.98,302.8,830.3,2356.,73.13,309.8,1034.,2532./
        DATA SA1C/.5,.5,.5,.5,7.571E-05,2.466E-04,8.113E-04,2.266E-03/
        DATA SA2C/.5,.5,.5,.5,5.508E-03,7.601E-03,7.880E-03,7.0096E-03/
C************************************************************************
        DATA IN_UNIT /20/                       ! INPUT UNIT NUMBER

        finter='data.req'
        print*,'Opening file data.req...'

 1110   if(finter .eq. ' ' .or. finter .eq. 'q') go to 1000
        open(10,file=finter,status='old',form='formatted',err=1111)
        go to  1112
 1111   print*,' '
        print*,'The file does not exist!'
 	print*,'Input another file containing the time intervals?'
        print*,' (usually data.req created after running request.exe)'
        read(*,'(a)') finter
        go to 1110

 1112   continue

c	open(unit=50,file='scratch:[schen]hkpwi.dat'
c	open(unit=50,file='hkpwi.dat'
c	1,form='unformatted',status = 'new')
	do i = 1,4
	   freqs(i) = saf1(i)
	   freqs(i+4) = saf2(i)
	   freqs(i+8) = sff(i)*1000
	   freqs(i+12) = sff(i+4)*1000
	   dfreqs(i) = sa1bw(i)
	   dfreqs(i+4) = sa2bw(i)
	   dfreqs(i+8) = sffbw(i)
	   dfreqs(i+12) = sffbw(i+4)
	enddo
c	write(50)freqs
c	write(50)dfreqs	   
c	write(*,*)freqs
c	write(*,*)dfreqs
c	print *,' type in time interval to average( units 5.75 sec)?'
c	read(*,*)iave
	iave=5.75
	tave = 5.7*(iave-1)
	dep = tave*1000d0
	bave = 0
	ibave = 0
C insert
C************************************************************************
C
C       OPEN THE HAWKEYE INPUT FILE             !/* VMS BLUE CODE */
C
C************************************************************************
500	JREC=0
	icount = 0
c	type *, 'Input filename''s yrdayhr'
c	read (*,*)idate
	read(10,*,end=1000) idate
	if( idate .lt. 0)goto 1000
        idespun =1
        do i=1,3
	   if(i.eq.1)write(filename,666)idate,'.despun'
	   if(i.eq.2)write(filename,666)idate,'.nondespun'
	   if(i.eq.2)write(filename,666)idate,'.dat'
           if(i.eq.2)idespun=0
           if(i.eq.3)idespun=2
           OPEN(UNIT=IN_UNIT,READONLY,STATUS='OLD',
	1       RECORDTYPE='VARIABLE',
	1       FORM='UNFORMATTED',file=filename,err=499)
	   print *,' file opened:',filename
           goto 600
499	   print *,' file not found: ',filename
	   if(i.eq.3)goto 500
        enddo
        
666	format('ndads::anon_dir:[data_dist.hawkeye]hk_',i7,a)
c666	format('scratch:[schen.ndads]hk_',i7,a)
C************************************************************************
c end insert

600 	continue

	write(filename2,667)idate,'.dat'
667	format('hkpwi_',i7,a)
	open(unit=50,file=filename2,form='unformatted',status = 'new')
	write(50)freqs
	write(50)dfreqs	   

   10   CONTINUE

C       READ IN ONE RECORD
        READ(IN_UNIT,END=500) dummy,ID
	if( cd(10) .le. 0.)then
	   print *,' cd(10)=',cd(10),'skip record'
	   goto 10
	endif
C
  20    CONTINUE    
        JREC=JREC+1
	IF( MOD(JREC,500) .EQ. 0)print *,' jrec=',jrec

	if( icount .eq. 0) then
C       TIME AND PASS NUMBER

        IYR=ID(3)                               !YEAR
        IDOY=ID(4)                              !DAY OF YEAR
        IHR=ID(5)                               !HOUR
        IMIN=ID(6)                              !MINUTE
        ISEC=ID(7)                              !SECOND
        IPASS=ID(11)                            !PASS NUMBER
	endif

	itod = (id(5)*60 + id(6))*6000 + id(7) !in centisec
	tod = itod/100.d0
	   ep = epoch0(id(3),id(4),tod)
	IF( ICOUNT .EQ. 0)THEN
	   krec = krec + 1
c	   type *, id(3),id(4),id(5),id(6),id(7),id(11)
	   itod0 = itod
	     ep0 = ep
c	   print *,' itod=',itod
	   tstart = tod	   
c	   print *,' krec, tstart=',krec, tstart
	   do i = 1,8
	   do j = 1,16
	      sa1(i,j) = 0.
	      sa2(i,j) = 0.
	      sfr(i,j) = 0.
	      sfr(i,j+16) = 0.
	    enddo
	    enddo
	endif
c	print *,' count=',icount,' tod=',tod
	call mag(bb0,bb1,jrec)
	if( bb0(3) .ne. -1.)then
	    bave = bave + (bb0(3) + bb0(6) + bb0(9))/3.
	    ibave = ibave + 1
	endif
c	print *,' bb0=',bb0
C************************************************************************
C       VLF DATA, RANGE, TEMPERATURE, MODE, AND CALIBRATION
C************************************************************************

 
        IR=1
        IF (ICOUNT.GT. 0) GO TO 301
        IVLF= (ID(12) .AND. 2**13 )/2**13       !ON/OFF
        IF (IVLF.EQ. 0) GO TO 10                !VLF NOT ON
        V(1)= (ID(12) .AND. 2**8 )/2**8	        !BAND WIDTH
        V(2)= (ID(12) .AND. 2**4 )/2**4         !MODE
        V(3)= (ID(12) .AND. 4)/4                !V DIFF GAIN
  301   CONTINUE
        J=MOD(ID(100),64)+1
        ICOUNT=ICOUNT+1
        IF (J.EQ.  9) V(4)=EW(2)                !TEMPERATURE
        IF (J.EQ.  9) IV(1)=IRAW(26)
        IF (J.EQ. 22) V(5)=EW(2)                !BOOM MOTOR TEMPERATURE
        IF (J.EQ. 22) IV(2)=IRAW(26)
        IF (J.EQ. 25) V(6)=EW(2)                !ANTENNA LENGTH WRONG
        IF (J.EQ. 25) IV(3)=IRAW(26)
        IF (J.EQ. 26) V(7)=EW(2)                !LOW VOLTAGE MONITOR
        IF (J.EQ. 26) IV(4)=IRAW(26)
        IF (J.EQ. 27) V(8)=EW(2)                !POWER SUPPLY STATUS
        IF (J.EQ. 27) IV(5)=IRAW(26)
        IF (J.EQ. 28) V(9)=EW(2)                !ELECTRIC ANTENNA LIMIT
        IF (J.EQ. 28) IV(6)=IRAW(26)
        IF (J.EQ. 29) V(10)=EW(2)               !TELESCOPIC BOOM LIMIT
        IF (J.EQ. 29) IV(7)=IRAW(26)
C
C       SPECTRUM ANALYZER DATA
C
        NJ=1+(ICOUNT-1)/4
        K=MOD(ID(100),4)
        IF (K.EQ. 0) L=2
        IF (K.EQ. 0) M=7
        IF (K.EQ. 1) L=8
        IF (K.EQ. 1) M=5
        IF (K.EQ. 2) L=6
        IF (K.EQ. 2) M=3
        IF (K.EQ. 3) L=4
        IF (K.EQ. 3) M=1
        if ( (k .lt. 0) .or. (k .gt. 3))then
            print *,' k out of range:',k
            print *,' id(100):',id(100)
        endif
        SA1(L,NJ)=IRAW(6)
        SA2(L,NJ)=IRAW(7)
        SA1(M,NJ)=IRAW(35)
        SA2(M,NJ)=IRAW(36)
C
C       STEP FREQUENCY RECEIVER ( SFR ) DATA
C
        MJ=1+(ICOUNT-1)/2
        K=MOD(ID(100),2)
        L=3
        M=5
        IF(K.EQ. 0) L=7
        IF(K.EQ. 0) M=1
        LL=L+1
        MM=M+1
        SFR(L,MJ)=IRAW(9)
        SFR(LL,MJ)=IRAW(8)
        SFR(M,MJ)=IRAW(38)
        SFR(MM,MJ)=IRAW(37)
C
C       ANTENNA POTENTIAL DIFFERENCE ( V-DIFF ) DATA
C
        IJ=ICOUNT
        VD(1,IJ)=IRAW(10)
        VD(2,IJ)=IRAW(22)
        VD(3,IJ)=IRAW(39)
        VD(4,IJ)=IRAW(43)
        VD(5,IJ)=IRAW(60)
        VD(6,IJ)=IRAW(64)
C

c        IF ((tod-tstart).LT. tave) GO TO 10             !TO OBTAIN ONE CALIBRATION
        IF ((ep-ep0).LT. dep) GO TO 10             !TO OBTAIN ONE CALIBRATION
c	print *,' raw sfr'
c	write(*,*) ((sfr(i,j),j=1,mj+1),i=1,8)	
C
C       CALCULATE ANTENNA LENGTH TIP TO TIP IN METERS
C
        X=IV(3)
        X=X*0.0196+0.005
        IF (X.GT. 1.4) GO TO 305
        XL=17.0736-1.1774*X
        XL=-9.7273+XL*X
        GO TO 306
  305   XL=16.1819-0.98509*X
        XL=-8.8570+XL*X
  306   CONTINUE
C
C       END LENGTH CALCULATION
C
C       CONVERT RAW DATA TO VOLTAGES
C
        AG=0.31342                              !GAIN OF AMPLIFIER
        VO=2.5497                               !OFF SET VOLTAGE
C
        IF (V(3).EQ. 1) AG=4.96819
        IF (V(3).EQ. 1) VO=2.5497
        DO 314 K=1,6
        VD0=VD0*0.0196+0.005
        IF (VD0.GT. 5.03) VD0=5.03
        VD0=-(VD0-VO)/AG                !CONVERT TO V rms
  314   CONTINUE
c        DO 315 J=1,64
c        DO 314 K=1,6
c        VD(K,J)=VD(K,J)*0.0196+0.005
c        IF (VD(K,J).GT. 5.03) VD(K,J)=5.03
c        VD(K,J)=-(VD(K,J)-VO)/AG                !CONVERT TO V rms
c  314   CONTINUE
c  315   CONTINUE
        DO 318 K=1,8
        DO 316 J=1,nj
        X=SA1(K,J)*0.0196+0.005
        IF (X.GT. 5.03) X=5.03
        CALL SSA1(X,K,Y)                        !CONVERT TO V rms
        SA1(K,J)=Y
        X=SA2(K,J)*0.0196+0.005
        IF (X.GT. 5.03) X=5.03
        CALL SSA2(X,K,Y)                        !CONVERT TO V rms
        SA2(K,J)=Y
  316   CONTINUE
        DO 317 J=1,mj
        X=SFR(K,J)*0.0196+0.005
        IF (X.GT. 5.03) X=5.03
        CALL SSFR(X,K,Y)                        !CONVERT TO V rms
	if( sfr(k,j) .eq. 0)then
c	   print *,' array el skipped (k,j)',k,j,' x,y=',x,y
	   y = -1
	endif
        SFR(K,J)=Y
  317   CONTINUE
  318   CONTINUE
C
c        DO 331 I=1,4                            !CONVERT VLF DATA TO
c        SA1C(I)=SA1C(I)*XL
c        SA2C(I)=SA2C(I)*XL                      !ELECTRIC AND MAGNETIC
c  331   CONTINUE
        SFRXL=XL/2.                             !SPECTRAL DENSITIES
        DO 333 II=1,8
        DO 332 JJ=1,16
        JK=JJ+16
        IF(II .LE. 4)THEN
        XX=SA1(II,JJ)/(SA1C(II)*xl)
	SA1(II,JJ)=XX*XX/SA1BW(II)
        XX=SA2(II,JJ)/(SA2C(II)*xl)
        SA2(II,JJ)=XX*XX/SA2BW(II)
	ELSE
        XX=SA1(II,JJ)/(SA1C(II))
	SA1(II,JJ)=XX*XX/SA1BW(II)
        XX=SA2(II,JJ)/(SA2C(II))
        SA2(II,JJ)=XX*XX/SA2BW(II)
	ENDIF
        XX=SFR(II,JJ)/SFRXL
        YY=SFR(II,JK)/SFRXL
	if( sfr(ii,jj) .eq. -1.)goto 332
        SFR(II,JJ)=XX*XX/SFFBW(II)
        SFR(II,JK)=YY*YY/SFFBW(II)
  332   CONTINUE
  333   CONTINUE
c average sfr, sa and output to file
        call average(sa1,sa2,sfr,spde,spdb)
	bave = bave/max(ibave,1)
	call trajectory(re,amlat,amlt,xgsm,ygsm,zgsm)  
c write output
	   tod = (ihr*60.d0 + imin)*60.d0 + isec/100.d0
	   epoch = epoch0(iyr,idoy,tod)
	   if( ifirst .eq. 0)then
	       ifirst = 1
	       write(50)iyr,idoy,ihr,imin,isec
	       write(*,*)iyr,idoy,ihr,imin,isec
	       write(*,*)' epoch=',epoch
	   endif
	WRITE(50)ipass,ep0 !YEAR,DAY,HRS,MIN,SEC*100
	WRITE(50)RE,AMLAT,AMLT,XGSM,YGSM,ZGSM 
	write(50)bave
	write(50)spde    !in units of v^2/m^2/Hz
	write(50)spdb    !in units of (nT)^2/Hz
        write(50)(id(12)-(id(12)/4)*4),mod(id(101),1024)
c reset parameter
	bave = 0
	ibave = 0
	icount = 0
	goto 10
C       END OF VLF DATA SECTION
1000	STOP
505     print *,' error encountered in file read, cd(10):',cd(10)
        goto 500
        END

	subroutine average(a,b,c,d,e)
	dimension a(8,16),b(8,16),c(8,32),d(16),e(8)
c average electric field spectral densities
	do i = 1,4
	   avea = 0
	   npta = 0
	   aveb = 0
	   nptb = 0
	do j = 1,16
	   if( a(i,j) .gt. 0.)then
	       avea = avea + a(i,j)
	       npta = npta + 1
	   endif
	   if( b(i,j) .gt. 0.)then
	       aveb = aveb + b(i,j)
	       nptb = nptb + 1
	   endif
	enddo
	   if( npta .ne. 0)avea = avea/npta
	   d(i) = avea
	   if( nptb .ne. 0)aveb = aveb/nptb
	   d(i+4) = aveb
	enddo
c sfr channels
	do i = 1,8
	   avec = 0
	   nptc = 0
	do j = 1,32
	   if( c(i,j) .gt. 0.)then
	       avec = avec + c(i,j)
	       nptc = nptc + 1
	   endif
	enddo
	   if( nptc .ne. 0)avec = avec/nptc
	   d(8+i) = avec
	enddo
c average magnetic field spectral densities
	do i = 5,8
	   avea = 0
	   npta = 0
	   aveb = 0
	   nptb = 0
	do j = 1,16
	   if( a(i,j) .gt. 0.)then
	       avea = avea + a(i,j)
	       npta = npta + 1
	   endif
	   if( b(i,j) .gt. 0.)then
	       aveb = aveb + b(i,j)
	       nptb = nptb + 1
	   endif
	enddo
	   if( npta .ne. 0)avea = avea/npta
	   e(i-4) = avea
	   if( nptb .ne. 0)aveb = aveb/nptb
	   e(i) = aveb
	enddo
	return
	end
C**********************************************************************
        SUBROUTINE SSA1(X,I,Y)
C**********************************************************************
C       ROUTINE TO CONVERT THE VLF LOG COMPRESSOR OUTPUT
C       TO THE CORRECTED INPUT VOLTAGES FOR THE SA-1 SYSTEM
C**********************************************************************
C       FREQUENCIES 1.78,5.62,17.8,56.2 HZ ELECTRIC AND MAGNETIC
C**********************************************************************
        DIMENSION B(8)
        DATA B/1.447,0.388,-0.006,-0.004,-0.133,-0.712,-0.746,-0.716/
        IF (X.LT. 4.601) GO TO 10
        Y=-5.03435+0.96335*X
        GO TO 20
 10     IF (X.LT. 4.192) GO TO 11
        Y=-6.32831+1.24456*X
        GO TO 20
 11     IF (X.LT. 3.599) GO TO 12
        Y=-4.18081+0.73229*X
        GO TO 20
 12     IF (X.LT. 3.081) GO TO 13
        Y=-5.62028+1.13227*X
        GO TO 20
 13     IF (X.LT. 2.524) GO TO 14
        Y=-4.33788+0.71605*X
        GO TO 20
 14     IF (X.LT. 1.988) GO TO 15
        Y=-5.35324+1.11837*X
        GO TO 20
 15     IF (X.LT. 1.401) GO TO 16
        Y=-4.57466+0.72671*X
        GO TO 20
 16     IF (X.LT. 1.004) GO TO 17
        Y=-5.30468+1.24766*X
        GO TO 20
 17     IF (X.LT. 0.414) GO TO 18
        Y=-4.88784+0.83264*X
        GO TO 20
 18     Y=-5.23333+1.66667*X
 20     Y=Y+B(I)
        Y=10.0**Y
        RETURN
        END
C**********************************************************************
        SUBROUTINE SSA2(X,I,Y)
C**********************************************************************
C       ROUTINE TO CONVERT THE VLF LOG COMPRESSOR OUTPUT
C       TO THE CORRECTED INPUT VOLTAGES FOR THE SA-2 SYSTEM
C**********************************************************************
C       FREQUENCIES 178,562,1.78K,5.62K HZ ELECTRIC AND MAGNETIC
C**********************************************************************
        DIMENSION B(8)
        DATA B/0.042,-0.004,-0.0055,0.093,-0.682,-0.75,-0.75,-0.6365/
        IF (X.LT. 4.641) GO TO 10
        Y=-5.27+X
        IF (I.EQ.4) Y=Y+0.083
        GO TO 20
 10     IF (X.LT. 4.155) GO TO 11
        Y=-6.12912+1.18510*X
        IF (I.EQ.4) Y=Y+0.048
        GO TO 20
 11     IF (X.LT. 3.623) GO TO 12
        Y=-4.09281+0.695018*X
        IF (I.EQ.4) Y=Y-0.022
        GO TO 20
 12     IF (X.LT. 3.087) GO TO 13
        Y=-5.67735+1.13235*X
        GO TO 20
 13     IF (X.LT. 2.527) GO TO 14
        Y=-4.36360+0.706838*X
        GO TO 20
 14     IF (X.LT. 2.007) GO TO 15
        Y=-5.44755+1.13578*X
        GO TO 20
 15     IF (X.LT. 1.391) GO TO 16
        Y=-4.60093+0.713992*X
        GO TO 20
 16     IF (X.LT. 0.943) GO TO 17
        Y=-5.32231+1.23250*X
        GO TO 20
 17     IF (X.LT. 0.486) GO TO 18
        Y=-4.84802+0.729694*X
        GO TO 20
 18     IF (X.LT. 0.264) GO TO 19
        Y=-5.03333+1.11111*X
        GO TO 20
 19     Y=-5.40+2.5*X
 20     Y=Y+B(I)
        Y=10.0**Y
        RETURN
        END
C************************************************************************
        SUBROUTINE SSFR(X,I,Y)
C************************************************************************
C       ROUTINE TO CONVERT THE VLF LOG COMPRESSOR OUTPUT
C       TO THE CORRECTED INPUT VOLTAGES FOR THE SFR SYSTEM
C************************************************************************
        DIMENSION B(8),A(8),C(8)
        DATA B/0.102,0.0188,-0.0069,-0.0315,0.0519,0.0169,0.0288,0.1383/
        DATA A/0.1452,0.0545,0.01815,0.0,0.0726,0.0363,0.0363,0.2087/
        DATA C/0.304,0.0967,0.0975,-0.0022,0.0553,-0.0029,0.0471,0.4526/
        IF (X.LT. 4.108) GO TO 10
        Y=-4.62046+0.90759*X
        Y=Y+A(I)
        IF (Y.LT. -0.345) GO TO 21
        Y=-11.95+2.5*X+C(I)
        GO TO 21
 10     IF (X.LT. 3.820) GO TO 11
        Y=-7.41225+1.58727*X
        GO TO 20
 11     IF (X.LT. 3.254) GO TO 12
        Y=-4.50012+0.82484*X
        GO TO 20
 12     IF (X.LT. 2.917) GO TO 13
        Y=-6.88878+1.55885*X
        GO TO 20
 13     IF (X.LT. 2.335) GO TO 14
        Y=-4.67732+0.80070*X
        GO TO 20
 14     IF (X.LT. 1.972) GO TO 15
        Y=-6.24591+1.47253*X
        GO TO 20
 15     IF (X.LT. 1.406) GO TO 16
        Y=-4.95012+0.81568*X
        GO TO 20
 16     IF (X.LT. 1.038) GO TO 17
        Y=-5.85304+1.45795*X
        GO TO 20
 17     IF (X.LT. 0.503) GO TO 18
        Y=-5.19802+0.82691*X
        GO TO 20
 18     Y=-7.7513E-06+4.82696E-05*X
c        Y=LOG10(Y)
	return
 20     Y=Y+B(I)
 21     Y=10.**Y
        RETURN
        END

	subroutine trajectory(re,amlat,amlt,xgsm,ygsm,zgsm) 

c	trj(4)  longitude in gse
c	trj(5)  latitude in gse
c	trj(6)  longitude in gmc converted to Magnetic local time
c       amlt =  mod( ( (trj(6)-trj(28)+180.)/15.+24.0 ),24.d0)
c	trj(7)  latitude in gmc
c	trj(9)  radius in RE
c	trj(10) gse-x
c	trj(11) gse-y
c	trj(12) gse-z
c	trj(13) gsm-x
c	trj(14) gsm-y
c	trj(15) gsm-z
        REAL*8 TRJ,CD
C**********************************************************************
        DIMENSION ID(324),TRJ(77),CD(21)
C**********************************************************************
        EQUIVALENCE (ID(130),CD(1))
        EQUIVALENCE (ID(180),TRJ(1))
	common/trecord/id
C
C*************************************************************************
C       TRAJECTORY INFORMATION
C*************************************************************************
C
	TRJ(9)=TRJ(9)/6378.16
	xgsm = trj(13)
 	ygsm = trj(14)
 	zgsm = trj(15)

    8   format(3x)
c  498   CONTINUE
c        WRITE (*,8)
c        WRITE (*,8)
c        WRITE (*,499)
c  499   FORMAT (30X,'TRAJECTORY DATA')
c  500   CONTINUE
c        WRITE (*,8)
c        WRITE (*,501)
c        WRITE (*,8)
c  501   FORMAT (27X,'TIME AND PASS NUMBER')
c        WRITE (*,502) ID(3),ID(4),ID(5),ID(6),ID(7),ID(11)
c        WRITE (*,8)
c  502   FORMAT (23X,6I6)
	amlt =  mod( ( (trj(6)-trj(28)+180.)/15.+24.0 ),24.d0)
	amlat = trj(7)  !gcc
	re = trj(9)

c        WRITE (*,503) TRJ(9),TRJ(4),TRJ(5)
c        WRITE (*,8)
c  503   FORMAT (23X,4HR = ,F7.3,3X,6HLON = ,F7.3,3X,6HLAT = ,F6.3)
c        WRITE (*,8)
c        WRITE (*,510)
c        WRITE (*,8)
c        WRITE (*,504) CD(10),CD(11),CD(12),CD(13),CD(14)
c  504   FORMAT (13X,5F10.3)
c        WRITE (*,8)
c        WRITE (*,505)
c  505   FORMAT (23X,'SATELLITE POSITION AND VELOCITY IN CI')
c        WRITE (*,8)
c        WRITE (*,506) TRJ(50),TRJ(51),TRJ(52),TRJ(53)
c  506   FORMAT(20X,'RAS =',F7.3,'DEC =',F7.3,'REV =',F7.3,'DEC =',F7.3)
c        WRITE (*,505)
c        WRITE (*,8)
c        WRITE (*,507)
c  507   FORMAT (35X,'GSE TO GSM')
c        WRITE (*,8)
c        WRITE (*,508) TRJ(32),TRJ(33),TRJ(34)
c        WRITE (*,508) TRJ(35),TRJ(36),TRJ(37)
c        WRITE (*,508) TRJ(38),TRJ(39),TRJ(40)
c  508   FORMAT (14X,3F15.6)
c        WRITE (*,8)
c        WRITE (*,509)
c  509   FORMAT (35X,'GEI TO GSE')
c        WRITE (*,8)
c        WRITE (*,508) TRJ(41),TRJ(42),TRJ(43)
c        WRITE (*,508) TRJ(44),TRJ(45),TRJ(46)
c        WRITE (*,508) TRJ(47),TRJ(48),TRJ(49)
c        WRITE (*,8)
c        WRITE (*,8)
c  510   FORMAT (31X,'ATTITUDE DATA')
	return
	end

	SUBROUTINE MAG(BCD2,BCD,IREC)
        INTEGER*4 ID,IRAW
C**********************************************************************
        REAL*8 TRJ,CD,EW,EM,EI,cd2(21)
C**********************************************************************
        DIMENSION ID(324),IRAW(64),CD(21),EW(2)
        DIMENSION BX(3),BY(3),BZ(3),BCD(9),bcd2(9)
        DIMENSION BB(3)
        DIMENSION EM(9),EI(9)
C**********************************************************************
        EQUIVALENCE (ID(34),IRAW(1)),(ID(104),EW(1)),(ID(130),CD(1))
        EQUIVALENCE (ID(242),EM(1)),(ID(260),EI(1))
	common/trecord/id
C************************************************************************
        IMAG = (ID(12) .AND. 2**16 )/2**16      !ON/OFF
c	PRINT *,' IMAG=',IMAG
        IF (IMAG.EQ. 0) return               !MAGNETOMETER NOT ON
        IX=ID(12)/4
        ICAL=0
        MAGRANGE =ID(12)-IX*4                   !RANGE 0 1 2 3
        JC=MOD(ID(101),1024)                    !S/C SEQUENCE NUMBER
	if( (jc .eq. 0) .or. (jc .eq. 1) ) then
	   do i = 1,9
	     bcd(i)=-1
	     bcd2(i)=-1
	   enddo
	   return
	endif	
        J=MOD (ID(100),32)+1                    !FRAME NUMBER
        IF (J.EQ. 6) SENTEMP=EW(2)
        IF (J.EQ. 7) ELETEMP=EW(2)
C
C       DESPIN CONVERTED MAGNETOMETER DATA
C
        XI=360.D0/CD(10)
        IF (CD(12).GE. 0.D0) GO TO 205          !FOR DESPUN DATA
        YI=irec-1                             !FOR NON-DESPUN DATA
        CD(12)=XI*YI*5.76D0                     !ASSUME START = 0.
 205    CONTINUE
       CD(1)=CD(1)+CD(12)+1.62D0*XI
       CD(1)=DMOD(CD(1),360.D0)
       CD(4)=CD(4)+CD(12)+3.51D0*XI
       CD(4)=DMOD(CD(4),360.D0)
       CD(7)=CD(7)+CD(12)+5.40D0*XI
       CD(7)=DMOD(CD(7),360.D0)
C
	do i = 1,21
	   cd2(i)=cd(i)
	enddo
       CALL ANGLE (EM,EI,CD)                   ! TRANSFORM TO GSM
C
C       PULL OUT RAW MAGNETOMETER DATA
C
        BX(1)=IRAW(19)
        BY(1)=IRAW(20)
        BZ(1)=IRAW(21)
        BX(2)=IRAW(40)
        BY(2)=IRAW(41)
        BZ(2)=IRAW(42)
        BX(3)=IRAW(61)
        BY(3)=IRAW(62)
        BZ(3)=IRAW(63)
c	print *,' mag raw =',bx(1),by(1),bz(1)
	btt = 257
	do i = 1,3
	if( bx(i) .ge. btt) goto 9898
	if( by(i) .ge. btt) goto 9898
	if( bz(i) .ge. btt) goto 9898
	enddo
	goto 9899
9898	   do i = 1,9
	     bcd(i)=-1
	     bcd2(i)=-1
	   enddo
	   return
9899	continue
        DO I=1,9
           BCD(I)=CD(I)  !3 measurements of B field in GSM coordinates
	ENDDO
 	DO K=1,3
	   BX(k)=BX(k)*0.0196+0.005            !CONVERT TO VOLTS
	   BY(k)=BY(k)*0.0196+0.005
	   BZ(k)=BZ(k)*0.0196+0.005
	ENDDO

        IF (MAGRANGE.EQ. 3) GO TO 260
        IF (MAGRANGE.EQ. 2) GO TO 250
        IF (MAGRANGE.EQ. 1) GO TO 240

	DO K=1,3
   	   BX(k)=BX(k)*59.996-149.962+3.685    !RANGE 1
	   BY(k)=BY(k)*60.200-150.438-0.908    !SCHONSTEDT
	   BZ(k)=BZ(k)*60.186-150.363          !CALIBRATIONS
	ENDDO
        GO TO 264
  240   CONTINUE
        DO K=1,3
	   BX(k)=BX(k)*180.565-451.345+3.685   !RANGE 2
	   BY(k)=BY(k)*179.854-449.488-0.908   !SCHONSTEDT
	   BZ(k)=BZ(k)*180.883-452.133         !CALIBRATIONS
	ENDDO
        GO TO 264
  250   CONTINUE
	DO K=1,3
	   BX(k)=BX(k)*597.133-1493.403+3.685  !RANGE 3
	   BY(k)=BY(k)*599.672-1498.909-0.908  !SCHONSTEDT
	   BZ(k)=BZ(k)*601.870-1503.170        !CALIBRATIONS
	ENDDO
        GO TO 264
  260   CONTINUE
	DO K=1,3
	   BX(k)=BX(k)*09977.3-24937.8+3.685   !RANGE 4
	   BY(k)=BY(k)*10003.7-24994.3-0.908   !SCHONSTEDT
	   BZ(k)=BZ(k)*10042.4-25100.9         !CALIBRATIONS
	ENDDO
  264   CONTINUE
	pi = 3.1415927
        arad = 180./pi
        DO K=1,3
	   brho =BX(k)*BX(k) + BY(k)*BY(k) 
	   X = brho + BZ(k)*BZ(k)
	   brho = sqrt(brho)
	   BB(k)=SQRT(X)
           theta = acos(bz(k)/bb(K))*arad
	   phi = acos(bx(k)/brho) 
	   if( by(k) .lt. 0) phi = 2*pi - phi
	   phi = phi*arad
	   cd2((k-1)*3 + 1) = phi
	   cd2((k-1)*3 + 2) = theta
	   cd2((k-1)*3 + 3) = bb(k)
	ENDDO
        CD2(1)=DMOD( CD2(1)+CD(12)+1.62D0*XI, 360.D0)
        CD2(4)=DMOD( CD2(4)+CD(12)+3.51D0*XI, 360.D0)
        CD2(7)=DMOD( CD2(7)+CD(12)+5.40D0*XI, 360.D0)
       CALL ANGLE (EM,EI,CD2)                   ! TRANSFORM TO GSM
        do k=1,9
	   bcd2(k)=cd2(k)
	enddo
c	print *,'     cd 			cd2'
c	do k=1,3
c	   print *, bcd((k-1)*3+1),bcd2((k-1)*3+1)
c	   print *, bcd((k-1)*3+2),bcd2((k-1)*3+2)
c	   print *, bcd((k-1)*3+3),bcd2((k-1)*3+3)
c	enddo
	RETURN
	END
