!======================================================================
!======================================================================
!				main lem_v5.f
!======================================================================
! lem_v5:
!	Galileo EPD Web Tool (LEMMS channels only)
!
! Features:
!	(1). Assumes that a user-generated, formatted, parameter 
!	     input file 'lem_v3' exists
!	(2). Assumes a Galileo NAV file with pre-determined contents
!	     See comments below.
!	(3). Assumes a EPD data file with pre-determined contents
!	(4). Restored (from lem_v1.f) useful diagnostic variables and
!	     pitch cosine in comoving frame
!	     See comments below.
!	(5). Code added to perform user-defined plasma convection
!	     velocity options jtvr=2
!
!	Fortran code written for Linux systems. Compile with:
!		g77 -o lem_v5 -O -W lem_v5.f
!	or
!		gfortran -o lem_v5 -O -W lem_v5.f
!
! History:
!	Created: 	23-Feb-07	rbd
!	Changed: 	06-Mar-07	rbd
!    v3 has TP, TO, and TS channels, 6/23/08, cp
!	Changed: 	11-Jul-08	rbd (see point (4) above)
!	Changed: 	24-Jul-08	rbd (flux -> rate_bgc, v4 ->v5)
!	Changed: 	28-Jul-08	rbd (see (5) above)
!======================================================================
!======================================================================
	program		lem_v5

	character*80	dfnav, dfdat, dfout

	integer*4  mdate(3), intim_nav(5), intim_dat(5), ms_rec(16)

	real*8	dsh_n, dsh_d, dsh_s, dsh_e

	real*4	tld(7),  ctld(7), stld(7)
        real*4	pld(16), cpld(16), spld(16)

	real*4	xld_mn(7,16), yld_mn(7,16), zld_mn(7,16)

	real*4	gal_jse(3), corot_jse(3), galv_sc(3) 
        real*4  corot_sc(3), jup_sc(3), jup_sc_unit(3), sun_sc(3)
	real*4	r_hat(3), c_hat(3), t_hat(3)

	real*4	pamag(16), bxd(16), byd(16), bzd(16)

! Auxilliary 16-element vector to hold input 16-sector ... 
!  ... background-corrected rate data (cnts/sec) and ...
	real*4	 A0_rate(16),A1_rate(16),A2_rate(16)
     	real*4	 A3_rate(16),A4_rate(16),A5_rate(16) 
	real*4   TP1_rate(16),TP2_rate(16),TP3_rate(16)
	real*4   TO2_rate(16),TO3_rate(16),TO4_rate(16)
	real*4   TS1_rate(16),TS2_rate(16),TS3_rate(16)
!  ... raw rate data (cnts/sec):
	real*4	 A0_rrate(16),A1_rrate(16),A2_rrate(16)
     	real*4	 A3_rrate(16),A4_rrate(16),A5_rrate(16) 
	real*4   TP1_rrate(16),TP2_rrate(16),TP3_rrate(16)
	real*4   TO2_rrate(16),TO3_rrate(16),TO4_rrate(16)
	real*4   TS1_rrate(16),TS2_rrate(16),TS3_rrate(16)

!======================================================================
! LEMMS:
!-----------------------
!	r_lem(i:15, k:16):
!		EPD LEMMS background-corrected rate channel i, sector k
!	rr_lem(i:15, j: 4, k:16):
!		EPD LEMMS raw rate channel i, species j, sector k
!----------------------------------------------------------------------
!       eg A0 is i=1, A1 is 2, A2 is 3, etc. 
!----------------------------------------------------------------------
	real*4	r_lem(15,16), rr_lem(15,16)	! 16-sect/15-ch (input)
	real*4	rav_lem(15)			! sect avg., 15 ch. (derived)
	real*4	rtf(15), egm(15), flx(15) 	! re-derived flux

	integer*4 nr_lem(15), nf_lem(15) 	! number good pnts. ch.
	real*4	  esp_ind(15)		 	! energy spectral indices

!----------------------------------------------------------------------
!  Geometric factors gf_lem(i:15) of LEMMS channel i:
!  Follow Mauk et al., JGR, 109, A09S12, 2004 and base TOF channels
!  on his formula in the appendix divided by 2.  If you divide formula
!  by 2.0, then formula and plot match.
!----------------------------------------------------------------------
	real*4	gf_lem(15)
	data	gf_lem/0.0052,0.0056,0.006, ! A0, A1, A2
     *		0.006,0.006,0.006,          ! A3, A4, A5
     *		0.0016,0.0014,0.0007,       ! TP1-TP3
     *		0.0028,0.0034,0.0037,       ! TO2-TO4
     *		0.0022,0.003,0.0035/        ! TS1-TS3

!----------------------------------------------------------------------
!  Energy passbands elo_lem(i:15, j:4) and ehi_lem(i:15, j:4) (keV) 
!  and efficiencies of LEMMS channel i, species j (2=H, 2=He, 3=O, 4=S)
!----------------------------------------------------------------------
	real*4	elo_lem(15,4), ehi_lem(15,4), eff_lem(15,4)
	real*4	atn(4), atm(4)
	data	atn/ 1.0, 2.0,  8.0, 16.0/		! atomic number
	data	atm/ 1.0, 4.0, 16.0, 32.0/		! mass number

!-----------------------
!  Species H (j=1):
!-----------------------
	data	elo_lem(1,1),	ehi_lem(1,1),	eff_lem(1,1)/	! A0/H
     *		22.0,		42.0,		1.00/,
     *		elo_lem(2,1),	ehi_lem(2,1),	eff_lem(2,1)/ 	! A1/H
     *		42.0,		65.0,		1.00/,
     *		elo_lem(3,1),	ehi_lem(3,1),	eff_lem(3,1)/	! A2/H
     *		65.0,		120.0,		1.00/,
     *		elo_lem(4,1),	ehi_lem(4,1),	eff_lem(4,1)/	! A3/H
     *		120.0,		280.0,		1.00/,
     *		elo_lem(5,1),	ehi_lem(5,1),	eff_lem(5,1)/	! A4/H
     *		280.0,		515.0,		1.00/,
     *		elo_lem(6,1),	ehi_lem(6,1),	eff_lem(6,1)/	! A5/H
     *		515.0,		825.0,		1.00/,
     *		elo_lem(7,1),	ehi_lem(7,1),	eff_lem(7,1)/	! TP1/H
     *		130.0, 		220.0,		1.0/
     *		elo_lem(8,1),	ehi_lem(8,1),	eff_lem(8,1)/	! TP2/H
     *		220.0, 		540.0,		1.0/
     *		elo_lem(9,1),	ehi_lem(9,1),	eff_lem(9,1)/	! TP3/H
     *		540.0, 		1140.0,		1.0/
     *		elo_lem(10,1),	ehi_lem(10,1),	eff_lem(10,1)/	! TO2/H
     *		0.0, 		0.0,		0.0/
     *		elo_lem(11,1),	ehi_lem(11,1),	eff_lem(11,1)/	! TO3/H
     *		0.0, 		0.0,		0.0/
     *		elo_lem(12,1),	ehi_lem(12,1),	eff_lem(12,1)/	! TO4/H
     *		0.0, 		0.0,		0.0/
     *		elo_lem(13,1),	ehi_lem(13,1),	eff_lem(13,1)/	! TS1/H
     *		0.0, 		0.0,		0.0/
     *		elo_lem(14,1),	ehi_lem(14,1),	eff_lem(14,1)/	! TS2/H
     *		0.0, 		0.0,		0.0/
     *		elo_lem(15,1),	ehi_lem(15,1),	eff_lem(15,1)/	! TS3/H
     *		0.0, 		0.0,		0.0/
!-----------------------
!-----------------------
!  Species He (j=2):
!-----------------------
	data elo_lem(1,2),	ehi_lem(1,2),	eff_lem(1,2)/ ! A0/He
     *		32.0,		50.0,		1.00/,
     *		elo_lem(2,2),	ehi_lem(2,2),	eff_lem(2,2)/ 	! A1/He
     *		50.0,		71.0,		1.00/,
     *		elo_lem(3,2),	ehi_lem(3,2),	eff_lem(3,2)/	! A2/He
     *		71.0,		129.0,		1.00/,
     *		elo_lem(4,2),	ehi_lem(4,2),	eff_lem(4,2)/	! A3/He
     *		129.0,		299.0,		1.00/,
     *		elo_lem(5,2),	ehi_lem(5,2),	eff_lem(5,2)/	! A4/He
     *		299.0,		560.0,		1.00/,
     *		elo_lem(6,2),	ehi_lem(6,2),	eff_lem(6,2)/	! A5/He
     *		560.0,		845.0,		1.00/,
     *		elo_lem(7,2),	ehi_lem(7,2),	eff_lem(7,2)/	! TP1/He
     *		0.0, 		0.0,		0.0/
     *		elo_lem(8,2),	ehi_lem(8,2),	eff_lem(8,2)/	! TP2/He
     *		0.0, 		0.0,		0.0/
     *		elo_lem(9,2),	ehi_lem(9,2),	eff_lem(9,2)/	! TP3/He
     *		0.0, 		0.0,		0.0/
     *		elo_lem(10,2),	ehi_lem(10,2),	eff_lem(10,2)/	! TO2/He
     *		0.0, 		0.0,		0.0/
     *		elo_lem(11,2),	ehi_lem(11,2),	eff_lem(11,2)/	! TO3/He
     *		0.0, 		0.0,		0.0/
     *		elo_lem(12,2),	ehi_lem(12,2),	eff_lem(12,2)/	! TO4/He
     *		0.0, 		0.0,		0.0/
     *		elo_lem(13,2),	ehi_lem(13,2),	eff_lem(13,2)/	! TS1/He
     *		0.0, 		0.0,		0.0/
     *		elo_lem(14,2),	ehi_lem(14,2),	eff_lem(14,2)/	! TS2/He
     *		0.0, 		0.0,		0.0/
     *		elo_lem(15,2),	ehi_lem(15,2),	eff_lem(15,2)/	! TS3/He
     *		0.0, 		0.0,		0.0/
!-----------------------
!-----------------------
!  Species O (j=3):
!-----------------------
	data	elo_lem(1,3),	ehi_lem(1,3),	eff_lem(1,3)/	! A0/O
     *		45.0,		107.0,		0.70/,
     *		elo_lem(2,3),	ehi_lem(2,3),	eff_lem(2,3)/ 	! A1/O
     *		93.0,		133.0,		0.80/,
     *		elo_lem(3,3),	ehi_lem(3,3),	eff_lem(3,3)/	! A2/O
     *		127.0,		210.0,		0.90/,
     *		elo_lem(4,3),	ehi_lem(4,3),	eff_lem(4,3)/	! A3/O
     *		200.0,		410.0,		1.00/,
     *		elo_lem(5,3),	ehi_lem(5,3),	eff_lem(5,3)/	! A4/O
     *		415.0,		650.0,		1.00/,
     *		elo_lem(6,3),	ehi_lem(6,3),	eff_lem(6,3)/	! A5/O
     *		650.0,		1000.0,		1.00/,
     *		elo_lem(7,3),	ehi_lem(7,3),	eff_lem(7,3)/	! TP1/O
     *		0.0, 		0.0,		0.0/
     *		elo_lem(8,3),	ehi_lem(8,3),	eff_lem(8,3)/	! TP2/O
     *		0.0, 		0.0,		0.0/
     *		elo_lem(9,3),	ehi_lem(9,3),	eff_lem(9,3)/	! TP3/O
     *		0.0, 		0.0,		0.0/
     *		elo_lem(10,3),	ehi_lem(10,3),	eff_lem(10,3)/	! TO2/O
     *		416.0, 		816.0,		1.0/
     *		elo_lem(11,3),	ehi_lem(11,3),	eff_lem(11,3)/	! TO3/O
     *		816.0, 		1792.0,		1.0/
     *		elo_lem(12,3),	ehi_lem(12,3),	eff_lem(12,3)/	! TO4/O
     *		1792.0, 	8992.0,		1.0/
     *		elo_lem(13,3),	ehi_lem(13,3),	eff_lem(13,3)/	! TS1/O
     *		0.0, 		0.0,		0.0/
     *		elo_lem(14,3),	ehi_lem(14,3),	eff_lem(14,3)/	! TS2/O
     *		0.0, 		0.0,		0.0/
     *		elo_lem(15,3),	ehi_lem(15,3),	eff_lem(15,3)/	! TS3/O
     *		0.0, 		0.0,		0.0/
!-----------------------
!-----------------------
!  Species S (j=4):
!-----------------------
	data	elo_lem(1,4),	ehi_lem(1,4),	eff_lem(1,4)/	! A0/S
     *		82.0,		157.0,		0.55/,
     *		elo_lem(2,4),	ehi_lem(2,4),	eff_lem(2,4)/ 	! A1/S
     *		137.0,		180.0,		0.72/,
     *		elo_lem(3,4),	ehi_lem(3,4),	eff_lem(3,4)/	! A2/S
     *		167.0,		270.0,		0.85/,
     *		elo_lem(4,4),	ehi_lem(4,4),	eff_lem(4,4)/	! A3/S
     *		250.0,		480.0,		1.00/,
     *		elo_lem(5,4),	ehi_lem(5,4),	eff_lem(5,4)/	! A4/S
     *		480.0,		750.0,		1.00/,
     *		elo_lem(6,4),	ehi_lem(6,4),	eff_lem(6,4)/	! A5/S
     *		750.0,		1150.0,		1.00/,
     *		elo_lem(7,4),	ehi_lem(7,4),	eff_lem(7,4)/	! TP1/S
     *		0.0, 		0.0,		0.0/
     *		elo_lem(8,4),	ehi_lem(8,4),	eff_lem(8,4)/	! TP2/S
     *		0.0, 		0.0,		0.0/
     *		elo_lem(9,4),	ehi_lem(9,4),	eff_lem(9,4)/	! TP3/S
     *		0.0, 		0.0,		0.0/
     *		elo_lem(10,4),	ehi_lem(10,4),	eff_lem(10,4)/	! TO2/S
     *		0.0, 		0.0,		0.0/
     *		elo_lem(11,4),	ehi_lem(11,4),	eff_lem(11,4)/	! TO3/S
     *		0.0, 		0.0,		0.0/
     *		elo_lem(12,4),	ehi_lem(12,4),	eff_lem(12,4)/	! TO4/S
     *		0.0, 		0.0,		0.0/
     *		elo_lem(13,4),	ehi_lem(13,4),	eff_lem(13,4)/	! TS1/S
     *		512.0, 		960.0,		1.0/
     *		elo_lem(14,4),	ehi_lem(14,4),	eff_lem(14,4)/	! TS2/S
     *		960.0, 		1984.0,		1.0/
     *		elo_lem(15,4),	ehi_lem(15,4),	eff_lem(15,4)/	! TS3/S
     *		1984.0, 	9920.0,		1.0/
!-----------------------
!======================================================================

	data	dtr, rtd/ 1.745329e-02, 5.729578e+01/
	data	pm, qe, sol/ 1.67e-24, 4.8e-10, 3.0e+10/
	data	c_ergtkev, c_kevterg/ 6.2422e+08, 1.602e-09/
	data	omega_J, rad_J/ 1.7585e-04, 7.14e+04/ 	! rad/sec, km
!
	npts1=0
	avg1=0.0
	sumsq1=0.0
	npts2=0.
	avg2=0.
	sumsq2=0.
!

!  Current date and time
	call idate (mdate)

!----------------------------------------------------------------------
! Open parameter input file (user-generated)
!----------------------------------------------------------------------
	open (unit=1, file='lem_v5.inp', form='formatted', status='old')

!----------------------------------------------------------------------
! Open Galileo navigation file:
!  nhnav = number of header info. records prior  to first nav record
!  dfnav = name of formatted nav file 
!----------------------------------------------------------------------
	do i=1,4
	   read (1,*)
	enddo
	read (1,*)  nhnav
	read (1,10) dfnav
  10	format (a80)
	open (unit=2, file=dfnav, form='formatted', status='old')
!  Read and discard nav file header records
	do i=1,nhnav
	  read (2,*)
	enddo

!----------------------------------------------------------------------
! Open EPD data file:
!  nhdat = number of header info. records prior to first data record
!  dfdat = name of formatted data file 
!----------------------------------------------------------------------
	read (1,*)
	read (1,*)
	read (1,*)  nhdat
	read (1,10) dfdat
	open (unit=4, file=dfdat, form='formatted', status='old')
!  Read and discard EPD data file header records
	do i=1,nhdat
	   read (4,*)
	enddo

!----------------------------------------------------------------------
! Enter name of output disk file for listing results:
!   ndo  = 0: do not write results to output file
!	   1: write results to output file
!  dfout = name of formatted output data file 
!----------------------------------------------------------------------
	read (1,*)
	read (1,*)
	read (1,*)  ndo
	read (1,10) dfout
	if (ndo .eq. 1) then
	   open (unit=8, file=dfout, form='formatted', status='new')
	endif

!----------------------------------------------------------------------
! Enter start and end times to process (all integers):
!  start time:	ksy, ksd, ksh, ksm, kss
!  end time:	key, ked, keh, kem, kes
!----------------------------------------------------------------------
	do i=1,3
	   read (1,*)
	enddo
	read (1,*) ksy, ksd, ksh, ksm, kss
	read (1,*)
	read (1,*) key, ked, keh, kem, kes

! Get decimal sequential hour form integer time values
	call c_seqt (1, dsh_s, idsh_s, ksy, ksd, ksh, ksm, kss)
	call c_seqt (1, dsh_e, idsh_e, key, ked, keh, kem, kes)

!----------------------------------------------------------------------
! Enter EPD channel (jcht) to tranform and ions species (jion):
!  jcht:  e.g., 1=A0, 2=A1, see input file for choices lem_v3.inp
!  jion:  1=H, 2=He, 3=O, 4=S
!----------------------------------------------------------------------
	do i=1,7
	   read (1,*)
	enddo
	read (1,*) jcht, jion

!----------------------------------------------------------------------
! Open file to write output from this run
!----------------------------------------------------------------------
!	open (unit=27, file='ratout.txt', form='formatted', status='new')
!	open (unit=28, file='gamout.txt', form='formatted', status='new')
!	open (unit=29, file='traout.txt', form='formatted', status='new')

	kk=0 ! initialize for write 

	
	if ((jcht .lt. 1) .or. (jcht .gt. 15)) then
	   print *, 'Invalid channel definition - stop.'
	   stop
	endif

	if ((jion .lt. 1) .or. (jion .gt. 4)) then
	   print *, 'Invalid species definition - stop.'
	   stop
	endif

! Check if channel 'jcht' has nonzero efficiency for species 'jion'
	if (eff_lem(jcht,jion) .lt. 0.10) then
	   print *, ' Species not measured by this channel - stop'
	   stop
	endif

!----------------------------------------------------------------------
! Get mean ion speed for this channel/species combination
!----------------------------------------------------------------------
	Tch = Sqrt(elo_lem(jcht,jion)*ehi_lem(jcht,jion))	! mean KE (keV)
	vch = Sqrt(2.0*Tch*c_kevterg/(atm(jion)*pm))/1.0e+05	! mean v, km/s

!	write(6,*) 'Channel number and kinetic energy:', jion, tch
!----------------------------------------------------------------------
! Energy Spectum Index Options:
!  d_sp_ind = default index to use when not calculable from EPD data
!----------------------------------------------------------------------
	do i=1,4
	   read (1,*)
	enddo
	read (1,*) d_esp_ind

!----------------------------------------------------------------------
! Plasma Convection Velocity Options (all speeds in km/s):
!  jvtr	: 0 = Rigid corotation direction and magnitude[*cfact=0.0-1.0]
!	: 1 = Rigid corotation direction, user-defined magnitude vc_cus
!	: 2 = User-defined convection velocity (vt_mag, vt_clt, vt_azm)
!
!  For option 2:
!	A local cartesian coordinate system centered on and moving with
!	Galileo is defined with:
!		c_hat[1:3] = unit vector in local corotation direction
!		r_hat[1:3] = unit radial vector from Jupiter to Galileo
!		t_hat[1:3] = unit vector c_hat x r_hat 
!	c_hat axis is polar axis from which polar angle vt_clt is defined
!		(vt_clt = [0,180] deg)
!	r_hat-t_hat plane is plane in which azimuth is measured, with 
!		vt_azm increasing from r_hat toward t_hat in right-hand 
!		sense along polar axis (vt_azm = [0,360))
!----------------------------------------------------------------------
	do i=1,6
	   read (1,*)
	enddo
	read (1,*) jvtr, cfact, vc_cus, vt_mag, vt_clt, vt_azm

	close (unit=1)
!
!	Now we have finished reading the input file

!----------------------------------------------------------------------
! Calculate look angles of EPD motor positions 1-7
! relative to spacecraft-frame cartesian system K[x,y,z]. 
!----------------------------------------------------------------------
	do m=1,7
	   tld(m)  = Float(m-1)*30.6
	   ctld(m) = Cos(tld(m)*dtr)
	   stld(m) = Sin(tld(m)*dtr)
	enddo

!----------------------------------------------------------------------
! Calculate look angles of EPD spin sectors 0-15
! relative to spacecraft-frame cartesian system K[x,y,z] 
!----------------------------------------------------------------------
	do n=1,16
	   pld(n)  = Float(2*n-1)*11.25		! 0-360 deg
	   cpld(n) = Cos(pld(n)*dtr)
	   spld(n) = Sin(pld(n)*dtr)
	enddo

!----------------------------------------------------------------------
! Calculate K[x,y,z] comps. of look-direction unit vector
! ld_mn = [xld_mn, yld_nm, zld_mn]. {Note: Unit vector of 
! incident particle is then -ld_mn.}
!----------------------------------------------------------------------
	do m=1,7
	   do n=1,16
	      xld_mn(m,n) = stld(m)*cpld(n)
	      yld_mn(m,n) = stld(m)*spld(n)
	      zld_mn(m,n) = ctld(m)
	   enddo
	enddo

	close (unit=1)			! close parameter input file

	n_nav = 0			! initialize nav record counter
	n_dat = 0			! initialize data record counter

!======================================================================
!			Nav Info. and EPD Data Input Loop
!======================================================================
 1000	continue

!----------------------------------------------------------------------
! Read a NAV file record (one record per time tag)
!----------------------------------------------------------------------
!	intim_nav(1) = year (4-digit) [i]
!	intim_nav(2) = doy (1-365 or 366) [i]
!	intim_nav(3) = hour (0-23) [i]
!	intim_nav(4) = min  (0-59) [i]
!	intim_nav(5) = sec  (0-59) [i]
!	fmpcn = current motor position (0-7) [r]
!	fmpln = motor position of previous spin (0-7) [r]
!	fnsbn = begin sector (0-63) [r]
!	fnsbn = end   sector (0-63) [r]
!	gal_jse[1:3]   = position vector of Galileo wrt Jupiter in 
!		         JSE system (km)
!	corot_jse[1:3] = corotation velocity at Galileo in JSE 
!			 system (km/s)
!	galv_sc[1:3]   = Galileo velocity in S/C system (km/s)
!	corot_sc[1:3]  = corotation unit vector direction at Galileo
!			 in S/C system
!	jup_sc[1:3]    = Position vector of Jupiter wrt Galileo in 
!		         S/C system (km)
!	sun_sc[1:3]    = Position vector of Sun wrt Galileo in 
!		         S/C system (km)
!----------------------------------------------------------------------
 1100 	read (2,*, end=2000) intim_nav, fmpcn, fmpln,
     *			     fnsbn, fnsen, gal_jse, corot_jse, galv_sc,
     *			     corot_sc, jup_sc, sun_sc
	mpcn = Ifix(fmpcn)		! current  motor position [i]
	mpln = Ifix(fmpln)		! previous motor position [i]
	nsbn = Ifix(fnsbn)		! start sector (0-63) [i]
	nsen = Ifix(fnsen)		! end   sector (0-63) [i]

	n_nav = n_nav + 1		! NAV record counter

	if (nsbn .ne. nsen) then
	   write (6,1110) n_nav
 1110	   format (' Start & end sectors unequal: NAV record ', i6/)
	   goto 1100
	endif

	kny = intim_nav(1)
	knd = intim_nav(2)
	knh = intim_nav(3)
	knm = intim_nav(4)
	kns = intim_nav(5)
	call c_seqt (1, dsh_n, idsh_n, kny, knd, knh, knm, kns)

!----------------------------------------------------------------------
! Read EPD data file records (multiple records per time tag)
!----------------------------------------------------------------------
!	intim_dat(1) = year [i]
!	intim_dat(2) = doy  [i]
!	intim_dat(3) = hour [i]
!	intim_dat(4) = min  [i]
!	intim_dat(5) = sec  [i]
!	mpcd = current motor positon (0-7) [i]
!	nsbd = begin sector of spin (0-63) [i]
!	nebd = end   sector of spin (0-63) [i]
!	pamag(16) = pitch angles of the 16 macro-sectors (0-180 deg) [r]
!	 *NOTE: The 16 pang elememts are always listed in order
!	        increasing macro-sector number 0-15
!	bxd(16) = x-component of the B (nT) wrt S/C system K[x,y,z]
!	byd(16) = y-component of the B (nT) wrt S/C system K[x,y,z]
!	bzd(16) = z-component of the B (nT) wrt S/C system K[x,y,z]
!----------------------------------------------------------------------
 1200 	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, pamag
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, bxd
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, byd
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, bzd
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A0_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A1_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A2_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A3_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A4_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A5_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TP1_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TP2_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TP3_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TO2_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TO3_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TO4_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TS1_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TS2_rate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TS3_rate

!
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  A0_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  A1_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  A2_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  A3_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  A4_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  A5_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  TP1_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  TP2_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  TP3_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  TO2_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  TO3_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  TO4_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  TS1_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  TS2_rrate
	read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed,  TS3_rrate
	

	kdy = intim_dat(1)
	kdd = intim_dat(2)
	kdh = intim_dat(3)
	kdm = intim_dat(4)
	kds = intim_dat(5)
	call c_seqt (1, dsh_d, idsh_d, kdy, kdd, kdh, kdm, kds)

	n_dat = n_dat + 1			! EPD data record counter

!----------------------------------------------------------------------
! Test for nav file and data file time synchonicity.
! Check for equivalence of integer year, doy, horr, min, and sec
! between the current NAV an EPD data record.
!----------------------------------------------------------------------

	if (dsh_n .lt. dsh_d) goto 1100
	if (dsh_n .gt. dsh_d) goto 1200

!----------------------------------------------------------------------
! Check nav file times against chosen start and end times.
!----------------------------------------------------------------------
	if (dsh_n .lt. dsh_s) goto 1000
	if (dsh_n .gt. dsh_e) goto 4000


!////////////////////////// SAMPLE OUTPUT /////////////////////////////
	write (6,1400)  intim_nav, n_nav,
     *			intim_dat, n_dat
 1400	format (/' nav: yrn,doyn,hhn,mmn,ssn,n_nav: ', 2i5,3i3,i5/
     *           ' dat: yrd,doyd,hhd,mmd,ssd,n_dat: ', 2i5,3i3,i5)

	gjr = Sqrt(gal_jse(1)**2 +   gal_jse(2)**2 +   gal_jse(3)**2)
	crj = Sqrt(corot_jse(1)**2 + corot_jse(2)**2 + corot_jse(3)**2)	
	crs = Sqrt(corot_sc(1)**2 +  corot_sc(2)**2 +  corot_sc(3)**2)
	gav = Sqrt(galv_sc(1)**2 +   galv_sc(2)**2 +   galv_sc(3)**2)
	jsc = Sqrt(jup_sc(1)**2 +    jup_sc(2)**2 +    jup_sc(3)**2)

! unit vector to Jupiter from Galileo = -unit radial from J to G 
	 jup_sc_unit(1) = jup_sc(1)/jsc
	 jup_sc_unit(2) = jup_sc(2)/jsc
	 jup_sc_unit(3) = jup_sc(3)/jsc

! test that unit vector from Galileo to Jupiter dotted with unit 
! corotation vector at Glileo vanishes (within single-precion)
	csc_d_jscu = 	corot_sc(1)*jup_sc_unit(1)
     *                + corot_sc(2)*jup_sc_unit(2)
     *                + corot_sc(3)*jup_sc_unit(3)
     
!----------------------------------------------------------------------
! rigid corotation speed in spin equatorial plane at radius gjr
	crv = omega_J*gjr*rad_J	
!----------------------------------------------------------------------

	write (6,1410) 	mpcn, mpln, nsbn, nsen, 
     *			gal_jse(1)*rad_J, 
     *			gal_jse(2)*rad_J,
     *			gal_jse(3)*rad_J,
     *			gjr*rad_J, galv_sc, gav, 
     *			corot_jse, crj, crv,
     *			corot_sc, crs, jup_sc_unit, 
     *                  csc_d_jscu, sun_sc,
     *			bxd, byd, bzd, A1_rate
 1410	format (' ', 4i5/
     *          '    gal_jse: ', 4(1pe11.3)/
     *          '    galv_sc: ', 4(1pe11.3)/
     *          '  corot_jse: ', 5(1pe11.3)/
     *          '   corot_sc: ', 4(1pe11.3)/
     *          'jup_sc_unit: ', 3(1pe11.3)/
     *          ' csc_d_jscu: ', 1(1pe13.5)/
     *          '     sun_sc: ', 3(1pe11.3)/
     *          '         Bx: ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '         By: ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '         Bz: ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '    A1_rate: ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3)/
     *          '             ', 4(1pe11.3))
!//////////////////////////////////////////////////////////////////////

!----------------------------------------------------------------------
! 			Process current record
!----------------------------------------------------------------------
!  Ignore data in motor position 0. Read in a new record.
	if (mpcn .eq. 0) goto 1000

!----------------------------------------------------------------------
! Fill in 2D arrays r_lem & rr_lem with vector (sectored) channel data
!----------------------------------------------------------------------
	do j=1,16
!   Sectored rates
	   r_lem(1,j)  =  A0_rate(j)
	   r_lem(2,j)  =  A1_rate(j)
	   r_lem(3,j)  =  A2_rate(j)
	   r_lem(4,j)  =  A3_rate(j)
	   r_lem(5,j)  =  A4_rate(j)
	   r_lem(6,j)  =  A5_rate(j)
	   r_lem(7,j)  =  TP1_rate(j)
	   r_lem(8,j)  =  TP2_rate(j)
	   r_lem(9,j)  =  TP3_rate(j)
	   r_lem(10,j) =  TO2_rate(j)
	   r_lem(11,j) =  TO3_rate(j)
	   r_lem(12,j) =  TO4_rate(j)
	   r_lem(13,j) =  TS1_rate(j)
	   r_lem(14,j) =  TS2_rate(j)
	   r_lem(15,j) =  TS3_rate(j)

!   Sectored fluxes
	   rr_lem(1,j)  =  A0_rrate(j)
	   rr_lem(2,j)  =  A1_rrate(j)
	   rr_lem(3,j)  =  A2_rrate(j)
	   rr_lem(4,j)  =  A3_rrate(j)
	   rr_lem(5,j)  =  A4_rrate(j)
	   rr_lem(6,j)  =  A5_rrate(j)
	   rr_lem(7,j)  =  TP1_rrate(j)
	   rr_lem(8,j)  =  TP2_rrate(j)
	   rr_lem(9,j)  =  TP3_rrate(j)
	   rr_lem(10,j) =  TO2_rrate(j)
	   rr_lem(11,j) =  TO3_rrate(j)
	   rr_lem(12,j) =  TO4_rrate(j)
	   rr_lem(13,j) =  TS1_rrate(j)
	   rr_lem(14,j) =  TS2_rrate(j)
	   rr_lem(15,j) =  TS3_rrate(j)
	enddo
!
!       Write the unprocessed intensity of the channel in question
!       to the output file number 27, named above
!       We do our own conversion from rate to intensity here
	fact1=ehi_lem(jcht,jion) - elo_lem(jcht,jion)
	fact2=gf_lem(jcht)*eff_lem(jcht,jion)*fact1
	do jj=1,16 
	if (r_lem(jcht,jj).gt.0) then
	aj=r_lem(jcht,jj)/fact2
!       This is the unprocessed intensity, aj
	write(27,*) n_dat,kdy,kdd,kdh,kdm,kds,aj
	npts1=npts1+1
!       we will average the logs for our metric so that
!       no point gets too much weight
	aja=alog10(aj)
	avg1=avg1+aja
	sumsq1=sumsq1+(aja*aja)
	endif
	enddo
!
!----------------------------------------------------------------------
!  Calculate sector-averaged rates and fluxes
!----------------------------------------------------------------------
	i1 = 1	! Have the capacity to do all channels
	i2 = 15 ! if wrong species will get NAN for tof channels

! Initialize accumulation buffers
	do i=i1,i2
	   rav_lem(i) = 0.0
	    nr_lem(i) = 0
	    nr_lem(i) = 0
	enddo

	do i=i1,i2
	   do j=1,16
!   Accumulate rates
	      rr = r_lem(i,j)
	      if ((rr .gt. 1.e-10) .and. (rr .lt. 1.e+10)) then
	         rav_lem(i) = rav_lem(i) + rr
	         nr_lem(i)  =  nr_lem(i) + 1
	      endif
	   enddo

!   Sector-averaged rate
	   if (nr_lem(i) .gt. 0) then
	      rav_lem(i) = rav_lem(i)/Float(nr_lem(i))
	   else
	      rav_lem(i) = -99.00
	   endif

!   Derive separate sector-averaged fluxes for chosen ion species (jion)
	   rtf(i) = gf_lem(i)*eff_lem(i,jion)*		! rate->flux
     *		    (ehi_lem(i,jion) - elo_lem(i,jion))	
	   egm(i) = Sqrt(ehi_lem(i,jion)*elo_lem(i,jion))	! log-mean

	   if (rav_lem(i) .gt. 0.0) then
	      flx(i) = rav_lem(i)/rtf(i)
	   else
	      flx(i) = -99.0
	   endif
	enddo

!======================================================================
!		Get estimates for energy spectral indices
!----------------------------------------------------------------------
!	Chans. A0-A5 (6 chans, 5 indices)
!		esp_ind(1) = Index using A0 & A1 [i=1 & i=2]
!		esp_ind(2) = Index using A1 & A2 [i=2 & i=3]
!		esp_ind(3) = Index using A2 & A3 [i=3 & i=4]
!		esp_ind(4) = Index using A3 & A4 [i=4 & i=5]
!		esp_ind(5) = Index using A4 & A5 [i=5 & i=6]
! temp		esp_ind(6) = defn.: esp_ind(5)
!----------------------------------------------------------------------
	do i=1,5
	   if ((flx(i) .gt. 0.0) .and. (flx(i+1) .gt. 0.0)) then
	      esp_ind(i) = -Log(flx(i)/flx(i+1))/Log(egm(i)/egm(i+1))
	   endif
	enddo
	esp_ind(6) = esp_ind(5)	! temp.: assign A5 the A4 spec. ind.
	
!----------------------------------------------------------------------
!	Chans. TP1-TP3 (3 chans., 2 indices)
!		esp_ind(7) = Index using TP1 & TP2 [i=7 & i=8]
!		esp_ind(8) = Index using TP2 & TP3 [i=8 & i=9]
! temp		esp_ind(9) = defn.: esp_ind(8)
!----------------------------------------------------------------------
	do i=7,8
	   if ((flx(i) .gt. 0.0) .and. (flx(i+1) .gt. 0.0)) then
	      esp_ind(i) = -Log(flx(i)/flx(i+1))/Log(egm(i)/egm(i+1))
	   endif
	enddo
	esp_ind(9) = esp_ind(8)	! temp.: assign TP3 the TP2 spec. ind.

!----------------------------------------------------------------------
!	Chans. TO2-TO4 (3 chans., 2 indices)
!		esp_ind(10) = Index using TO2 & TO3 [i=10 & i=11]
!		esp_ind(11) = Index using TO3 & TO4 [i=11 & i=12]
! temp		esp_ind(12) = defn.: esp_ind(11)
!----------------------------------------------------------------------
	do i=10,11
	   if ((flx(i) .gt. 0.0) .and. (flx(i+1) .gt. 0.0)) then
	      esp_ind(i) = -Log(flx(i)/flx(i+1))/Log(egm(i)/egm(i+1))
	   endif
	enddo
	esp_ind(12) = esp_ind(11) ! temp.: assign TO4 the TO3 spec. ind.

!----------------------------------------------------------------------
!	Chans. TS1-TS3 (3 chans., 2 indices)
!		esp_ind(13) = Index using TS1 & TS2 [i=13 & i=14]
!		esp_ind(14) = Index using TS2 & TS3 [i=14 & i=15]
! temp		esp_ind(15) = defn.: esp_ind(14)
!----------------------------------------------------------------------
	do i=13,14
	   if ((flx(i) .gt. 0.0) .and. (flx(i+1) .gt. 0.0)) then
	      esp_ind(i) = -Log(flx(i)/flx(i+1))/Log(egm(i)/egm(i+1))
	   endif
	enddo
	esp_ind(15) = esp_ind(14) ! temp.: assign TS3 the TS2 spec. ind.

!       write gamma values used to file 28, named above
!	if not species in question will get NAN for egms
!
!	print *, '    esp_ind: ', (esp_ind(i), i=i1,i2)
!======================================================================


!======================================================================
!		Frame Transformation Procedure
!======================================================================

!----------------------------------------------------------------------
! Determine Plasma Convection Velocity
!----------------------------------------------------------------------
	Vcon_x = 0.0
	Vcon_y = 0.0
	Vcon_z = 0.0
!----------------------------------------------------------------------
!  jvtr: 0 = Rigid corotation direction and magnitude[*cfact=0.0-1.0]
!----------------------------------------------------------------------
	if (jvtr .eq. 0) then
	   Vcon_x  = corot_sc(1)*crv*cfact - galv_sc(1)
	   Vcon_y  = corot_sc(2)*crv*cfact - galv_sc(2)
	   Vcon_z  = corot_sc(3)*crv*cfact - galv_sc(3)

!----------------------------------------------------------------------
!  jvtr: 1 = Rigid corotation direction, user-defined magnitude vc_cus
!----------------------------------------------------------------------
	else if (jvtr .eq. 1) then
	   Vcon_x = corot_sc(1)*vc_cus - galv_sc(1)
	   Vcon_y = corot_sc(2)*vc_cus - galv_sc(2)
	   Vcon_z = corot_sc(3)*vc_cus - galv_sc(3)

!----------------------------------------------------------------------
!  jvtr: 2 = User-defined convection velocity (vt_mag, vt_clt, vt_azm)
!----------------------------------------------------------------------
!  Construct local cartensian system centered on and moving with S/C:
!	r_hat[1:3] = unit radial vector from Jupiter to Galileo
!	c_hat[1:3] = unit vector in direction of corotation
!	t_hat[1:3] = unit vector c_hat x r_hat
! 	Let c_hat axis be polar axis from which polar angle (colatitude)
!	   angle vt_clt [0,180] is defined.
!	Let R_hat-t_hat plane be equatorial plane in which azimuth angle
!	   vt_azm [0,360) is defined, with vt_azm measured from r_hat 
!	   CCW toward t_hat [note: wrt usual spherical polars, identify
!	   z -> c-hat, x -> r_hat, y -> t_hat]
!----------------------------------------------------------------------
	else if (jvtr .eq. 2) then
	   do i=1,3
	      c_hat(i) =  corot_sc(i)		! unit corotation vector
	      r_hat(i) = -jup_sc_unit(i)	! unit radial vector J->G
	   enddo

! construct unit vector t_hat (=c_hat x r_hat) normal to c_hat-r_hat plane
	   t_hat(1) = c_hat(2)*r_hat(3) - c_hat(3)*r_hat(2)
	   t_hat(2) = c_hat(3)*r_hat(1) - c_hat(1)*r_hat(3) 
	   t_hat(3) = c_hat(1)*r_hat(2) - c_hat(2)*r_hat(1) 

! get components of convection velocity along local axes [c_hat,r_hat,t_hat]
	   V_c_hat = vt_mag*Cos(vt_clt*dtr)
	   V_r_hat = vt_mag*Sin(vt_clt*dtr)*Cos(vt_azm*dtr)
	   V_t_hat = vt_mag*Sin(vt_clt*dtr)*Sin(vt_azm*dtr)

! get components of convection velocity along S/C axes [xsc,ysc,zsc]
	   Vcx = V_c_hat*c_hat(1) + V_r_hat*r_hat(1) + V_t_hat*t_hat(1)
	   Vcy = V_c_hat*c_hat(2) + V_r_hat*r_hat(2) + V_t_hat*t_hat(2)
	   Vcz = V_c_hat*c_hat(3) + V_r_hat*r_hat(3) + V_t_hat*t_hat(3)

! subtract off S/C velocity
	   Vcon_x = Vcx - galv_sc(1)
	   Vcon_y = Vcy - galv_sc(2)
	   Vcon_z = Vcz - galv_sc(3)
	endif

	Vcon    = Sqrt(Vcon_x**2 + Vcon_y**2 + Vcon_z**2)
	Vcon_ux = Vcon_x/Vcon		! x-comp. of unit vector Vcon_u
	Vcon_uy = Vcon_y/Vcon		! y-comp. of unit vector Vcon_u
	Vcon_uz = Vcon_z/Vcon		! z-comp. of unit vector Vcon_u

	print *, 'Vcon_x,Vcon_y,Vcon_z,Vcon: ', 
     *            Vcon_x,Vcon_y,Vcon_z,Vcon
!
!       Write a header for 28, gamma file
	if (kk.eq.0) write(28,*) jcht, jion, crv, vcon
	kk=1
	write(28,*) n_dat,kdy,kdd,kdh,kdm,kds,esp_ind(jcht)
!----------------------------------------------------------------------
!  Get begin sector ns_beg (0-15) from micro-sector nsbd (0-63) by
!  doing integer divide by 4 (i.e., 4 micro-sectors per macro-sector)
!----------------------------------------------------------------------
	ns_beg = nsbd/4				! start sector (0-15)
	print *, 'nsbd, ns_beg:', nsbd, ns_beg

!----------------------------------------------------------------------
!  Now calculate the 16 consecutive sectors in current record
!----------------------------------------------------------------------
	do n=1,16
	   mtmp = ns_beg + (n-1)
	   ms_rec(n) = Mod(mtmp,16)		! must be [1,16]
	   print *, 'n, ms_rec(n):', n, ms_rec(n)
	enddo

!----------------------------------------------------------------------
!  For current record, determine look directions of the 16 sector 
!  positions at current motor position. Look direction for current 
!  record, sector n [1-16], is unit vector ld = (xld,xld,xld), and 
!  vector velocity of entering ion is v = (vxld,vyld,vzld).
!----------------------------------------------------------------------

!  Current motor positon [1-7]
	mp = mpcn
	print *, 'mp,tld(mp):', mp, tld(mp)

!
!----------------------------------------------------------------------
!  Look directions 16 sectors of current record	
!  ms_rec(n:16) contains numbers of the 16 macro-sectors [0,15]
!----------------------------------------------------------------------
	do n=1,16
	   ns = ms_rec(n) + 1		! index [1-16] of sector [0-15]

!----------------------------------------------------------------------
!  EPD look direction unit vector: 
!   LEMMS in 180 deg end: multiply components by -1
!----------------------------------------------------------------------
	   xld = -xld_mn(mp,ns)	! x-comp. wrt SC sys. K[x,y,z]
	   yld = -yld_mn(mp,ns)	! y-comp. wrt SC sys. K[x,y,z]
	   zld = -zld_mn(mp,ns)	! z-comp. wrt SC sys. K[x,y,z]

!----------------------------------------------------------------------
!  Unit vector of particle entry velocity in SC frame K[x,y,z]:
!		t_v_mn = [t_vx_mn, t_vy_mn, t_vz_mn]
!----------------------------------------------------------------------
	   t_vx_mn = -xld	! vx-comp. wrt SC sys. K[x,y,z]
	   t_vy_mn = -yld	! vy-comp. wrt SC sys. K[x,y,z]
	   t_vz_mn = -zld	! vz-comp. wrt SC sys. K[x,y,z]

!----------------------------------------------------------------------
!  Cos(t_v_mn,Vcon_u): Cos[angle between (unit vector of particle 
!  entry velocity in SC frame K) and (unit vector of convection 
!  velocity Vcon_u)]
!----------------------------------------------------------------------
	   Cos_tv_Vcon  = t_vx_mn*Vcon_ux + 
     *			  t_vy_mn*Vcon_uy + 
     *			  t_vz_mn*Vcon_uz
	   Sin2_tv_Vcon = 1.0 - Cos_tv_Vcon**2	   

!----------------------------------------------------------------------
!  Derived particle speed in SC frame K[x,y,z] at this direction (m,n) 
!  given the look-direction-independent constant speed vch defined
!  in comoving frame K'[x',y',z']:
!----------------------------------------------------------------------
	   v_mn = vch*((Vcon/vch)*Cos_tv_Vcon +
     *                  Sqrt(1.0 - ((Vcon/vch)**2)*Sin2_tv_Vcon))
c	   v_mn = vch + Vcon*Cos_tv_Vcon	! approx. to O(Vcon/vch)

!----------------------------------------------------------------------
!  Derived components of unit vector to look-direction-independent 
!  constant speed vch defined in comoving frame K'[x',y',z']:
!		t_vp_mn = [t_vxp_mn, t_vyp_mn, t_vzp_mn]
!  NOTE: Needed to calc. pitch cosine in comoving frame K'[x',y',z']
!----------------------------------------------------------------------
	   t_vxp_mn = (t_vx_mn*v_mn - Vcon_x)/vch
	   t_vyp_mn = (t_vy_mn*v_mn - Vcon_y)/vch
	   t_vzp_mn = (t_vz_mn*v_mn - Vcon_z)/vch

!----------------------------------------------------------------------
!  Dot product of ld=(xld,yld,zld) with B=(bxd,byd,bzd) gives the 
!  pitch cosine of the sector's look direction (sld). To get the pitch
!  cosine of a particle that enters that sector, change the sign
!  of the sector's look direction pitch cosine.
!----------------------------------------------------------------------
	   Bx   = bxd(n)
	   By   = byd(n)
	   Bz   = bzd(n)
	   Bm   = Sqrt(Bx**2 + By**2 + Bz**2)
	   b_ux = Bx/Bm
	   b_uy = By/Bm
	   b_uz = Bz/Bm
	   
!----------------------------------------------------------------------
!  Pitch cosine and pitch angle of particle velocity 
!  in SC frame K[x,y,z]:
!----------------------------------------------------------------------
	   pcos_v = t_vx_mn*b_ux + t_vy_mn*b_uy + t_vz_mn*b_uz
	   pang_v = ACos(pcos_v)*rtd

!----------------------------------------------------------------------
!  Pitch cosine and pitch angle of particle velocity
!  in comoving frame K'[x',y',z']:
!----------------------------------------------------------------------
	   pcos_vp = t_vxp_mn*b_ux + t_vyp_mn*b_uy + t_vzp_mn*b_uz
	   pang_vp = ACos(pcos_vp)*rtd

!...................................................................
! Note: Because the NAV file always lists the magnetic pitch-angles
!  	of the 16 macro-sectors in order, from sector 0 through 15,
!	to compare my results pang to those on the NAV records use
!	the index ns to increment pamag, i.e., as pamag(ns).
!...................................................................
c	   print *, '    n,ms_rec(n),ns,pld(ns),pang_v,pamag(ns):',
c     *                  n,ms_rec(n),ns,pld(ns),pang_v,pamag(ns)

	   rate = r_lem(jcht,n)			! sector rate from file
	   flxd = rr_lem(jcht,n)			! sector flux from file 
	   flux = rate/rtf(jcht)		! rederived flux

!----------------------------------------------------------------------
!  Velocity components & K.E. of incident particle in SC frame K
!----------------------------------------------------------------------
	   vx_mn = -t_vx_mn*vch		! x-comp.
	   vy_mn = -t_vy_mn*vch		! y-comp.
	   vz_mn = -t_vz_mn*vch		! z-comp.
	   T_mn  = vx_mn**2 + vy_mn**2 + vz_mn**2
	   v_mn  = Sqrt(T_mn)

!----------------------------------------------------------------------
!  Velocity components & K.E. of incident particle in comoving frame K'
!----------------------------------------------------------------------
	   vxp_mn = vx_mn - Vcon_x	! x'-comp
	   vyp_mn = vy_mn - Vcon_y	! y'-comp
	   vzp_mn = vz_mn - Vcon_z	! z'-comp
	   Tp_mn  = vxp_mn**2 + vyp_mn**2 + vzp_mn**2	   
	   vp_mn  = Sqrt(Tp_mn)

!----------------------------------------------------------------------
!  Get energy spetral index for this channel.
!   Try indices eps_ind(jcht) and esp_ind(jcht+1) first;
!   use default d_esp_ind if these are unsatisfactory.
!   Note: ==> For now, channels A1-A5 only, i1=2, i2=6
!----------------------------------------------------------------------
	   sp_ind = 0.0  ! for the compiler
	   if ((jcht .ge. i1) .and. (jcht .lt. i2)) then
	      if (Abs(esp_ind(jcht)) .lt. 6.0) then
	         sp_ind = esp_ind(jcht)
c	      elseif (Abs(esp_ind(jcht+1)) .lt. 6.0) then
c	         sp_ind = esp_ind(jcht+1)
	      else
	         sp_ind = d_esp_ind
	      endif
	   endif
	  
	   cgf    = (Tp_mn/T_mn)**(sp_ind + 1) 	! C-G factor
	   rate_p = rate*cgf
	   flux_p = flux*cgf

!	   print *,
!     *     'n,ms_rec(n),pamag(ns),pang_v,v_mn,vp_mn,',
!     *     'sp_ind,rate,rate_p:',
!     *      n,ms_rec(n),pamag(ns),pang_v,v_mn,vp_mn,
!     *      sp_ind,rate,rate_p
c       This is what Rob had:
c	   if (ndo .gt. 0) then
c	      write (8,*) intim_dat, n, ms_rec(n), pamag(ns), pang_v,
c     *		          v_mn, vp_mn, rate, rate_p
c	   endif
	   if (ndo .gt. 0) then
	      write (8,*) intim_dat(2),intim_dat(3),intim_dat(4),
     *     intim_dat(5),pang_v, pang_vp, rate, rate_p
	   endif

!
!       write the transformed intensity to file 29, named above
!       If the gamma is 1-3, there is little reason
!       to preserve pitch angle and do the spectral slopes separately
!
	if (flux_p.gt.0.) then
	write(29,*) n_dat,kdy,kdd,kdh,kdm,kds,flux_p
	npts2=npts2+1
!       Average the log of the new transformed intensity
	fpa=alog10(flux_p)
	avg2=avg2+fpa
	sumsq2=sumsq2+(fpa*fpa)
	endif
!
	enddo



!======================================================================

	goto 1000

 2000	write (6, 2010)
 2010	format (' EOF on Nav data file.')
	goto 5000

 3000	write (6, 3010)
 3010	format (' EOF on EPD data file.')
	goto 5000

 4000	write (6, 4010)
 4010	format (' User-defined end time encountered.')

 5000	close (unit = 2)
 	close (unit = 4)

	if (npts1.gt.0) then
	avg1=avg1/float(npts1)
	fact3=sumsq1-(npts1*avg1*avg1)
	sdev1=sqrt(fact3/float(npts1-1))
	write(6,*) 0.,npts1,avg1,sdev1,sdev1/avg1
	endif
	if (npts2.gt.0) then
	avg2=avg2/float(npts2)
	fact4=sumsq2-(npts2*avg2*avg2)
	sdev2=sqrt(fact4/float(npts2-1))
	write(6,*) 1.,npts2,avg2,sdev2,sdev2/avg2
	endif

	write(6,*) 'jcht, jion, crv, vcon'
	write(6,*) jcht, jion, crv, vcon

	stop
	end


!====================================================================
! subroutine c_seqt: 
!--------------------
! Sequential time in hours or days is elapsed time from
! 1960/001/00:00:00
!
!	k=-1:	Convert decimal sequential HOUR (seqt) to integer:
!		sequential HOUR (iseqt), year (iyr), doy (idoy),
!		hour (ihr), minute (imn), and second (isc).
!
!	k=-2:	Convert decimal sequential DAY (seqt) to integer:
!		sequential DAY (iseqt), year (iyr), doy (idoy),
!		hour (ihr), minute (imn), and second (isc).
!
!	k= 1:	Convert integer year (iyr), doy (idoy), hour (ihr),
!		minute (imn), and second (isc) to integer sequential
!		HOUR (iseqt) and decimal sequential HOUR (seqt).
!
!	k= 2:	Convert integer year (iyr), doy (idoy), hour (ihr),
!		minute (imn), and second (isc) to integer sequential
!		DAY (iseqt) and decimal sequential DAY (seqt).
!
! notes:
!	(1) uses 4-digit year, iyr (1960 .ge. iyr .le. 2022)
!	(2) idoy = 1 to 365 or 366
!	(3) ihr  = 0 to 23
!	(4) imn  = 0 to 59
!	(5) isc  = 0 to 59
!	(6) seqt must be real*8 in calling program
!
! created: 04-jun-97	RBDecker
! changed: 28-jun-00	RBDecker
!====================================================================

	subroutine	c_seqt (k,seqt,iseqt,iyr,idoy,ihr,imn,isc)
	integer*4	jsd(64)
	real*8		seqt, sqd, f, fhr, fmn, fsc, eps

! doy 0, 1960 through 2022 (add 1 to get doy 1)
	data jsd/    0,   366,   731,  1096,	! 000/1960 to 000/1963
     *		  1461,  1827,  2192,  2557,	! 000/1964 to 000/1967
     *		  2922,  3288,  3653,  4018,	! 000/1967 to 000/1971
     *		  4383,  4749,  5114,  5479,	! 000/1972 to 000/1975
     *		  5844,  6210,  6575,  6940,	! 000/1976 to 000/1979
     *		  7305,  7671,  8036,  8401,	! 000/1980 to 000/1983
     *		  8766,  9132,  9497,  9862,	! 000/1984 to 000/1987
     *		 10227, 10593, 10958, 11323,	! 000/1988 to 000/1991
     *		 11688, 12054, 12419, 12784,	! 000/1992 to 000/1995
     *		 13149, 13515, 13880, 14245,	! 000/1996 to 000/1999
     *		 14610, 14976, 15341, 15706,	! 000/2000 to 000/2003
     *		 16071, 16437, 16802, 17167,	! 000/2004 to 000/2007
     *		 17532, 17898, 18263, 18628,	! 000/2008 to 000/2011
     *		 18993, 19359, 19724, 20089,	! 000/2012 to 000/2015
     *		 20454, 20820, 21185, 21550,	! 000/2016 to 000/2019
     * 		 21915, 22281, 22646, 23011/	! 000/2020 to 000/2023
	data eps/ 1.0d-04/

! Compiler wants next two statements
	sqd  = 0.0d0
	isqd = 0

!--------------------------------------------------------------------
! k=-1 or -2: seqt -> [iseqt=Int(seqt), iyr, idoy, ihr, imn, isc]
!--------------------------------------------------------------------
	if (k .lt. 0) then
	   iseqt = seqt				! integer sequential time

	   if (k .eq. -1) then
	      sqd  = seqt/24.d0			! decimal sequential day, k=-1
	      isqd = sqd			! integer sequential day, k=-1
	   else if (k .eq. -2) then
	      sqd  = seqt			! decimal sequential day, k=-2
	      isqd = sqd			! integer sequential day, k=-2
	   end if

	   nn = 1
	   do n=1,63
	      if ((isqd .gt. jsd(n)) .and. (isqd .le. jsd(n+1))) then
	         nn = n
	         go to 10
	      end if
	   end do

 10	   iyr  = nn + 1959
	   idoy = isqd - jsd(nn)
	   fhr  = (sqd - Float(isqd))*24.d0
	   ihr  =  fhr + eps			! integer hour of doy
	   fmn  = (fhr - Float(ihr))*60.d0
	   imn  =  fmn + eps			! integer minute of hour
	   fsc  = (fmn - Float(imn))*60.d0	
	   isc  =  fsc  + eps			! integer second of minute

! test for isc=60, imn=60, ihr=24
	   if (isc .lt. 60) go to 20
	   isc  = 0
	   imn  = imn + 1
 20	   if (imn .lt. 60) go to 25
	   imn  = 0
	   ihr  = ihr + 1
 25	   if (ihr .lt. 24) return
	   ihr  = 0
	   isqd = isqd + 1
	end if

!--------------------------------------------------------------------
! k=1 or 2: [iyr, idoy, ihr, imn, isc] -> [seqt, iseqt=Int(seqt)]
!--------------------------------------------------------------------
	if (k .gt. 0) then
	   jyr = iyr
	   if (jyr .lt. 100) jyr = jyr + 1900
	   n = jyr - 1959			! no. years since 1960

	   isqd = jsd(n) + idoy			! integer sequential doy
	   f    =   Float(ihr)
     *		  + Float(imn)/60.d0
     *		  + Float(isc)/3600.d0		! decimal hour of day

	   if (k .eq. 1) then
	      khr   = isqd*24
	      iseqt = khr + ihr			! integer sequential hour
	      seqt  = Float(khr) + f		! decimal sequential hour
	   else if (k .eq. 2) then
	      iseqt = isqd			! integer sequential day
	      seqt  = Float(iseqt) + f/24.d0	! decimal sequential day
	   end if

	end if

	return
	end

