c............ Turbulence_Hourly_Analysis.for ..........................
c
c  Routine to read merged MAG/SWEPAM data and generate hourly
c    averages of basic solar wind statistics.
c
c  C.W.S. -- 11/10/2005 & 5/30/2006.
c.........................................................................................................................
c
c	Wind Attributes		     Variable    Represented By
c   						
c
c	Solar Wind Velocity          Vp(60)      data(i,4)
c	Solar Wind Velocity (r)      Vr(60)      data(i,1)
c	Solar Wind Velocity (t)      Vt(60)      data(i,2)
c	Solar Wind Velocity (n)      Vn(60)      data(i,3)
c	Solar Wind Velocity (x)      Vx(60)      
c	Solar Wind Velocity (y)      Vy(60)
c	Solar Wind Velocity (z)      Vz(60)
c	Proton Temperature           Tp(60)      data(i,11)
c	Proton Density	             Np(60)      data(i,10)
c	Magnetic Field (r)           Br(60)      data(i,5)
c	Magnetic Field (t)           Bt(60)      data(i,6)
c       Magnetic Field (n)           Bn(60)      data(i,7)
c	Magnetic Field (x)           Bx(60)      
c	Magnetic Field (y)           By(60)  
c	Magnetic Field (z)           Bz(60)
c	RMS Magnetic Field           Brms(60)    data(i,9)
c	Magnetic Field Magnitude     Bmag(60)    data(i,8)
c	Alfven Speed  	             Va(60)      data(i,15)     
c       Anisotropic Magnetic Field   Baniso
c	Anisotropic Velocity	     Vaniso
c	plasma beta	             beta  	 data(i,12)
c	alpha			     alpha       data(i,13)
c				     zlambda     data(i,14)
c			
c	    
c	Statistical Attributes			  Variable
c
c       <dB^2>                                    Bvar
c	Uncertainty in <db^2>                     Bvar_sigma
c	<dV^2>	                                  Vvar
c	Uncertainty in <dV^2>			  Vvar_sigma
c	Correlation between V and B fluctuations  VdotB    
c	delta_Br Squared		          delta_Br2
c	delta_Bt Squared			  delta_Bt2
c	delta_Bn Squared			  delta_Bn2
c	delta_Vr Squared			  delta_Vr2
c	delta_Vt Squared			  delta_Vt2
c	delta_Vn Squared			  delta_Vn2
c	V dot B					  VdotB
c	V dot B average				  VdotB_aver
c	V dot B instantaneous			  VdotB_inst
c	V dot B average uncertainty		  VdotB_aver_var
c	V dot B instantaneous uncertainty         VdotB_inst_var			  
c
c	var2					  standard deviation
c	var4                                      variance                                     
c				    
c	Time Attributes	
c
c	iyear_start = Start year of data analysis
c	iyear_stop = End year of data analysis
c	time_start = Start time of data anlysis [days]
c	time_stop = Stop time of data analysis [days]
c
c
c	Other Attributes		Variable
c	
c	Spacecraft Location (x)		Xsc(60)
c	Spacecraft Location (y)		Ysc(60)
c	Spacecraft Location (z)		Zsc(60)
c	
c.........................................................................................................................
c
c.. native data declarations:                        
	integer*4 icnt_VdotB                                                   
	real*8 VdotB, VdotB2, VdotB_var, VdotB_temp
	real*8 Baniso, Vaniso, EVi, EBi, Btrace(60), Vtrace(60) 
	real*8 S2_Bx(60), S2_By(60), S2_Bz(60),S2_Vx(60), S2_Vy(60), S2_Vz(60), S2_Np(60), S2_Tp(60)
	real*8 S2_Tp2(60),S2_Np2(60),S2_Bx2(60),S2_By2(60),S2_Bz2(60),S2_Vx2(60),S2_Vy2(60),S2_Vz2(60)
	real*8 sigma_NpStruc(60),sigma_TpStruc(60),sigma_BxStruc(60),sigma_ByStruc(60),sigma_BzStruc(60)
	real*8 sigma_VxStruc(60),sigma_VyStruc(60),sigma_VzStruc(60),sigma_BtraceStruc(60),sigma_VtraceStruc(60)
	real*8 data(60,15), y(60), ymean, yvar2, yvar4, mean(60)
	real*8 var2(60), var4(60), delta(60,15), delta2(60,15),delta_dat(60)
	real*8 var2_mnfld(60), var4_mnfld(60), mean_mnfld(60)
	real*8 one_hour
c
c.. data declarations for data input:
	integer*4 npts,nptss
	integer*4 iyear, icnt_mag(60)
	real*8 time,times, Np(60), Tp(60), alpha(60), Vp(60), Vr(60)
	real*8 Vt(60), Vn(60)
	real*8 Br(60), Bt(60), Bn(60), Bmag(60), zlambda(60), Brms(60)
	real*8 delta_Vr(60), delta_Vt(60), delta_Vn(60), delta_Br(60)
	real*8 delta_Bt(60), delta_Bn(60)
	real*8 Bx(60), By(60), Bz(60), Vx(60), Vy(60), Vz(60)
	real*8 Xsc(60), Ysc(60), Zsc(60)
c
c...Means and variance from the MNFLD routine
	real*8 Bx_mean,Bx_mean2,Bx_mean4,By_mean,By_mean2,By_mean4,Bz_mean,Bz_mean2,Bz_mean4
     	real*8 Vx_mean,Vx_mean2,Vx_mean4,Vy_mean,Vy_mean2,Vy_mean4,Vz_mean,Vz_mean2,Vz_mean4
c

	common /SWEPAM/ Bx, By, Bz, Vx, Vy, Vz
	common /hourly/ time, Np, Tp, alpha, Vp, Vr, Vt, Vn,
     .		Br, Bt, Bn, Bmag, zlambda, delta_dat, Xsc, Ysc, Zsc, 
     .		icnt_mag,times,iyear
	common /structure/ S2_Tp, S2_Np, S2_Bx, S2_By, S2_Bz, S2_Vx,S2_Vy, S2_Vz,k,N,Btrace,Vtrace,
     .                     S2_Tp2, S2_Np2, S2_Bx2, S2_By2, S2_Bz2, S2_Vx2, S2_Vy2, S2_Vz2,sigma_VxStruc,
     .                     sigma_VyStruc,sigma_VzStruc,sigma_NpStruc,sigma_TpStruc,sigma_VtraceStruc,
     .                     sigma_BxStruc,sigma_ByStruc,sigma_BzStruc,sigma_BtraceStruc
	common /mnfld_data/ Bx_mean,Bx_mean2,Bx_mean4,By_mean,By_mean2,By_mean4,Bz_mean,Bz_mean2,Bz_mean4,
     .			    Vx_mean,Vx_mean2,Vx_mean4,Vy_mean,Vy_mean2,Vy_mean4,Vz_mean,Vz_mean2,Vz_mean4
c
c.. data declarations for data selection:
	integer*4 iyear_start, iyear_stop
	real*8 time_start, time_stop
	character*40 merged_files_listing
	common /timeselect/ iyear_start, iyear_stop, time_start,
     .                      time_stop,merged_files_listing
c
dd	real*8
	real*8 Tp_A,Tp_sigmaA,Tp_B,Tp_sigmaB,Np_A,Np_sigmaA,Np_B
	real*8 Np_sigmaB,Bx_A,Bx_sigmaA,Bx_B,Bx_sigmaB,By_A,By_sigmaA
	real*8 By_B,By_sigmaB,Bz_A,Bz_sigmaA,Bz_B,Bz_sigmaB,Vx_A
	real*8 Vx_sigmaA,Vx_B,Vx_sigmaB,Vy_A,Vy_sigmaA,Vy_B,Vy_sigmaB
	real*8 Vz_A,Vz_sigmaA,Vz_B,Vz_sigmaB,Btrace_A, Btrace_B
	real*8 Btrace_sigmaA,Btrace_sigmaB
     	real*8 Vtrace_A, Vtrace_B,Vtrace_sigmaA,Vtrace_sigmaB
	real*8 Bx_A_rtn,Bx_sigmaA_rtn,Bx_B_rtn,Bx_sigmaB_rtn,By_A_rtn,By_sigmaA_rtn
	real*8 By_B_rtn,By_sigmaB_rtn,Bz_A_rtn,Bz_sigmaA_rtn,Bz_B_rtn,Bz_sigmaB_rtn,Vx_A_rtn
	real*8 Vx_sigmaA_rtn,Vx_B_rtn,Vx_sigmaB_rtn,Vy_A_rtn,Vy_sigmaA_rtn,Vy_B_rtn,Vy_sigmaB_rtn
	real*8 Vz_A_rtn,Vz_sigmaA_rtn,Vz_B_rtn,Vz_sigmaB_rtn,Btrace_A_rtn,Btrace_B_rtn
	real*8 Btrace_sigmaA_rtn,Btrace_sigmaB_rtn
     	real*8 Vtrace_A_rtn,Vtrace_B_rtn,Vtrace_sigmaA_rtn,Vtrace_sigmaB_rtn
c
c...Passfit exchanges data from the linefit
	common /passfit/ Tp_A,Tp_sigmaA,Tp_B,Tp_sigmaB,Np_A,Np_sigmaA,
     .			 Np_B,Np_sigmaB,
     .                   Bx_A,Bx_sigmaA,Bx_B,Bx_sigmaB,By_A,By_sigmaA,
     .                   By_B,By_sigmaB,
     .                   Bz_A,Bz_sigmaA,Bz_B,Bz_sigmaB,Vx_A,Vx_sigmaA,
     .                   Vx_B,Vx_sigmaB,
     .                   Vy_A,Vy_sigmaA,Vy_B,Vy_sigmaB,Vz_A,Vz_sigmaA,
     .                   Vz_B,Vz_sigmaB,
     .                   Btrace_A,Btrace_sigmaA,Btrace_B,
     .                   Btrace_sigmaB,
     .                   Vtrace_A,Vtrace_sigmaA,Vtrace_B,Vtrace_sigmaB
c
c.. data declarations for housekeeping variables:
	integer*4 iy1, iy2,i,j,k,N       
	real*8 bad, pi, t1, t2
	common /house/ bad, pi, t1, t2, iy1, iy2
c
c...data declarations for detrend
	integer*8 idetrend
c
	bad = -9999.9
	pi = 4.*atan(1.)
	one_hour = 1./24.
c
	Do 5 i=1,60
	   Bx(i) = bad
	   By(i) = bad
	   Bz(i) = bad
	   Vx(i) = bad
	   Vy(i) = bad
	   Vz(i) = bad
	   Br(i) = bad
	   Bn(i) = bad
	   Bt(i) = bad
 5	continue

c
c
c....open linefit files
c      open(unit=90,file='/home/students/tessein/windstats/output/linefit_Tp.out',form='formatted',status='new')
c      open(unit=91,file='/home/students/tessein/windstats/output/linefit_Np.out',form='formatted',status='new')
c      open(unit=92,file='/home/students/tessein/windstats/output/linefit_Bx.out',form='formatted',status='new')
c      open(unit=93,file='/home/students/tessein/windstats/output/linefit_By.out',form='formatted',status='new')
c      open(unit=94,file='/home/students/tessein/windstats/output/linefit_Bz.out',form='formatted',status='new')
c      open(unit=95,file='/home/students/tessein/windstats/output/linefit_Vx.out',form='formatted',status='new')
c      open(unit=96,file='/home/students/tessein/windstats/output/linefit_Vy.out',form='formatted',status='new')
c      open(unit=97,file='/home/students/tessein/windstats/output/linefit_Vz.out',form='formatted',status='new')
c      open(unit=98,file='/home/students/tessein/windstats/output/linefit_Btr.out',form='formatted',status='new')
c      open(unit=99,file='/home/students/tessein/windstats/output/linefit_Vtr.out',form='formatted',status='new')
c       	
c.. start it up:
	write(6,90)
   90	format(/' Routine to read merged MAG/SWEPAM data files ',
     .	       /' and analyze hourly turbulence parameters.')
c
	write(6,95)
   95	format(' Give start and stop times for analysis: (iyr, day)')
	read(5,*) iyear_start, time_start, iyear_stop, time_stop
	write(6,*) iyear_start, time_start, iyear_stop, time_stop
c
c.. get file listing:
	write(6,100)
  100	format(' Give local filename with dsnames listing merged MAG/SWEPAM data (a40):')
	read(5,110) merged_files_listing
  110	format(a40)
	write(6,120) merged_files_listing
  120	format(' File name is ',a40)
c
  130	idetrend = 0
	write(6,140)
  140	format(' Do you want to detrend? (1 = yes)')
	read(5,*) idetrend
	if(idetrend .eq. 0) write(6,*)' Choice is Do not Detrend'
	if(idetrend .eq. 1) write(6,*)' Choice is Detrend'
	if((idetrend .gt. 1) .or. (idetrend .lt. 0)) then
		write(6,*) 'Try Again'
		goto 130
	endif
c.. open output file:
	open(unit=18, file='THA.out', form='formatted', status='old')
c
c.. open file listing and get first data file:
	open(unit=14, file=merged_files_listing, form='formatted', status='old')
	call getfile(iflag,iwhy)
	if(iflag .eq. -1) then
	    write(6,*) '  Did not find any useful data in specified time range - STOP!'
	    close(unit=14)
	    stop
	endif
c...establish time-keeping rules
	iy1 = iyear_start
	t1 = time_start
	t2 = time_start + one_hour
	iy2 = iyear_start
	if(((iy1/4)*4 .ne. iy1) .and. (t2 .gt. 366.)) then
	    iy2 = iy2 + 1
	    t2 = t2 - 366.
	endif
	if(((iy1/4)*4 .eq. iy1) .and. (t2 .gt. 367.)) then
	    iy2 = iy2 + 1
	    t2 = t2 - 367.
	endif
c
  200	continue				! loop point for hourly analysis	
	Do 201 i=1,60
	   Bx(i) = bad
	   By(i) = bad
	   Bz(i) = bad
	   Vx(i) = bad
	   Vy(i) = bad
	   Vz(i) = bad
	   Br(i) = bad
	   Bn(i) = bad
	   Bt(i) = bad
	   Vr(i) = bad
	   Vt(i) = bad
	   Vn(i) = bad
 201	continue


c.. get an hour of merged data:
	call gethour(iflag,iwhy,npts,Brms)
c
	if(iflag .lt. 0) go to 1000
c.. compute statistics:			! this is where the work is done!
c  
c....Bring the data into one array, making it more managable
	do 210 i=1,npts
	  data(i,1) = Vr(i)
	  data(i,2) = Vt(i)
	  data(i,3) = Vn(i)
	  data(i,4) = Vp(i)
	  data(i,5) = Br(i)
	  data(i,6) = Bt(i)
	  data(i,7) = Bn(i)
	  data(i,8) = Bmag(i)
	  data(i,9) = Brms(i) 
	  data(i,10) = Np(i)
	  data(i,11) = Tp(i)
	  data(i,12) = alpha(i)
	  if((abs(Np(i)-bad) .gt. 1.) .and. (abs(Tp(i)-bad) .gt. 1.) 
     .      .and. (abs(Bmag(i)-bad) .gt. 1.)) then		! beta
	      data(i,13) = 4.03E-1 * Np(i)*Tp(i)*0.86E-4/Bmag(i)**2
	   else
              data(i,13) = bad
	  endif
	  if(abs(Np(i)-bad) .gt. 1.) then			! zlambda_ii
	      data(i,14) = 228. / sqrt(Np(i))
	    else
		data(i,14) = bad
	  endif
	  if((abs(Np(i)-bad).gt.1.).and.(abs(Bmag(i)-bad).gt.1.))then
	      data(i,15) = 21.8 * ( Bmag(i) / sqrt(Np(i)))	! Va
	    else
		data(i,15) = bad
	  endif
  210	continue
c
	nptss = npts
	do 220 j=1,15
	  do 215 i=1,npts
	    y(i) = data(i,j)
  215	  continue
	  call meanvar(y, nptss, ymean, yvar2, yvar4)
	  mean(j) = ymean
	  var2(j) = yvar2
	  var4(j) = yvar4
  220	continue
c
	Bvar = var2(5)**2 + var2(6)**2 + var2(7)**2
	sigma_Bvar = sqrt(var4(5)**2 + var4(6)**2 + var4(7)**2)
	Vvar = var2(1)**2 + var2(2)**2 + var2(3)**2
	sigma_Vvar = sqrt(var4(1)**2 + var4(2)**2 + var4(3)**2)
c
	do 450 j=1,15
	  if(abs(mean(j)-bad) .gt. 1.) then
		do 400 i=1,npts
		  if(abs(data(i,j)-bad) .gt. 1.) then
			delta(i,j) = data(i,j) - mean(j)
			delta2(i,j) = delta(i,j)**2
		    else
			delta(i,j) = bad
			delta2(i,j) = bad
		  endif
  400		continue
	    else
	        do 440 i=1,npts
		    delta(i,j) = bad
		    delta2(i,j) = bad
  440	        continue
	  endif
  450	continue
c
	icnt_VdotB = 0
	VdotB_aver = 0.
	VdotB_aver2 = 0.
	VdotB_inst = 0.
	VdotB_inst2 = 0.
	do 470 i=1,npts 	                     
	if((abs(delta(i,1)-bad) .gt. 1.) .and. (abs(delta(i,5)-bad) .gt. 1.) .and.
     .     (abs(delta(i,2)-bad) .gt. 1.) .and. (abs(delta(i,6)-bad) .gt. 1.) .and.
     .     (abs(delta(i,3)-bad) .gt. 1.) .and. (abs(delta(i,7)-bad) .gt. 1.) .and.
     .	   (abs(delta(i,3)-bad) .gt. 1.) .and. (abs(delta(i,7)-bad) .gt. 1.)) then
	   VdotB_temp=(delta(i,1)*delta(i,5)+delta(i,2)*delta(i,6)+delta(i,3)*delta(i,7)) 
		icnt_VdotB = icnt_VdotB + 1
		VdotB_aver = VdotB_aver + VdotB_temp
		VdotB_aver2 = VdotB_aver2 + VdotB_temp**2
		EVi = delta(i,1)**2 + delta(i,2)**2 + delta(i,3)**2
		EBi = delta(i,5)**2 + delta(i,6)**2 + delta(i,7)**2
		VdotB_temp = VdotB_temp / (sqrt(EVi)*sqrt(EBi))
		VdotB_inst = VdotB_inst + VdotB_temp
		VdotB_inst2 = VdotB_inst2 + VdotB_temp**2
	  endif
  470	continue
	if(icnt_VdotB .gt. 1) then
	    zcnt_VdotB = float(icnt_VdotB)
	    VdotB_aver = VdotB_aver / zcnt_VdotB
	    VdotB_aver2 = VdotB_aver2 / zcnt_VdotB
	    VdotB_aver_var = sqrt(VdotB_aver2 - VdotB_aver**2)
	    VdotB_aver = VdotB_aver / (sqrt(Vvar)*sqrt(Bvar))
	    VdotB_aver2 = VdotB_aver2 / (sqrt(Vvar)*sqrt(Bvar))
	    VdotB_aver_var = VdotB_aver_var / (sqrt(Vvar)*sqrt(Bvar))
	    VdotB_inst = VdotB_inst / zcnt_VdotB
	    VdotB_inst2 = VdotB_inst2 / zcnt_VdotB
	    VdotB_inst_var = sqrt(VdotB_inst2 - VdotB_inst**2)
	  else
	    VdotB_aver = bad
	    VdotB_aver2 = bad
	    VdotB_aver_var = bad
	    VdotB_inst = bad
	    VdotB_inst2 = bad
	    VdotB_inst_var = bad
	endif
c
c.....Compute anisotropic B and v
	Do 500 i=1,npts
	 Bx(i) = Br(i)
	 By(i) = Bt(i)
	 Bz(i) = Bn(i)
	 Vx(i) = Vr(i)
	 Vy(i) = Vt(i)
	 Vz(i) = Vn(i)
 500	continue
	if(idetrend .eq. 1) call detrend(npts)
	call structurefunc(npts) !compute structure functions for RTN coordinates
c.....compute Baniso and Vaniso
c
c...Put results of RTN linefit into different variables before being overwritten
	Btrace_A_rtn = Btrace_A 
	Btrace_B_rtn = Btrace_B
	Btrace_sigmaA_rtn = Btrace_sigmaA
	Btrace_sigmaB_rtn = Btrace_sigmaB
        Vtrace_A_rtn = Vtrace_A
	Vtrace_B_rtn = Vtrace_B
	Vtrace_sigmaA_rtn = Vtrace_sigmaA
	Vtrace_sigmaB_rtn = Vtrace_sigmaB
	Bx_A_rtn = Bx_A
	Bx_B_rtn = Bx_B
	Bx_sigmaA_rtn = Bx_sigmaA
	Bx_sigmaB_rtn = Bx_sigmaB
        By_A_rtn = By_A
        By_B_rtn = By_B
        By_sigmaA_rtn = By_sigmaA
	By_sigmaB_rtn = By_sigmaB
	Bz_A_rtn = Bz_A
	Bz_B_rtn = Bz_B
	Bz_sigmaA_rtn = Bz_sigmaA
        Bz_sigmaB_rtn = Bz_sigmaB
	Vx_A_rtn = Vx_A
	Vx_B_rtn = Vx_B
	Vx_sigmaA_rtn = Vx_sigmaA
	Vx_sigmaB_rtn = Vx_sigmaB
        Vy_A_rtn = Vy_A
        Vy_B_rtn = Vy_B
	Vy_sigmaA_rtn = Vy_sigmaA
	Vy_sigmaB_rtn = Vy_sigmaB
	Vz_A_rtn = Vz_A
	Vz_B_rtn = Vz_B
	Vz_sigmaA_rtn = Vz_sigmaA
	Vz_sigmaB_rtn = Vz_sigmaB          
c
	if(idetrend .eq. 1) then
	  Do 501 i=1,npts
	   Bx(i) = Br(i)
	   By(i) = Bt(i)
	   Bz(i) = Bn(i)
	   Vx(i) = Vr(i)
	   Vy(i) = Vt(i)
	   Vz(i) = Vn(i)
 501	  continue
	endif
	call mnfld(npts,Baniso,Vaniso) !rotate coordinate system to obtain xyz coordinates 
	if(idetrend .eq. 1) call detrend(npts)
c
c....compute structure functions on xyz coordinates
	call structurefunc(npts)    
c
c
c
c
c
c ---------- lots more to do here as we evolve the project !!!!
c
c
c.. output:
c
c...first block
c...Vp mean, Vp variance, Vr mean, Vr variance
c...Vt mean, Vt variance, Vn mean, Vn variance
c...Va mean, Va variance
c...Tp mean, Tp variance, Np mean,Np variance
c...beta mean, beta variance, zlambda_ii mean, zlambda_ii variance
c...Bmag mean, Bmag variance, Brms mean, Brms variance
c...Br mean, Br variance, Bt mean, Bt variance
c...Bn mean, Bn variance, Bvar, sigma_Bvar
c...Bmag variance**2, Bmag**2 variqnce, Br variance**2, Br**2 variance
c...Bt variance**2, Bt**2 variance, Bn variance**2, Bn**2 variance
c...Vp variance**2, Vp**2 variance, Vr variance**2, Vr**2 variance
c...Vt variance**2, Vt**2 variance, Vn variance**2, Vn**2 variance

c
c...second block
c...V variance,sigma V variance  
c...VdotB average,VdotB average variance, VdotB instantaneous,VdotBinstantaneous variance  
c...B anisotropy,V anisotropy  
c...Btrace A,Btrace B,Btrace A sigma,Btrace B sigma ! Linefit data for y = a + bx
c...Vtrace A,Vtrace B,Vtrace A sigma,Vtrace B sigma
c...Temp A,Temp A sigma,Temp B,Temp B sigma
c...Np A,Np A sigma,Np B,Np B sigma
c...Bx A,Bx A sigma,Bx B,Bx B sigma
c...By A,By A sigma,By B,By B sigma
c...Bz A,Bz A sigma,Bz B,Bz B sigma
c...Vx A,Vx A sigma,Vx B,Vx B sigma
c...Vy A,Vy A sigma,Vy B,Vy B sigma
c...Vz A,Vz A sigma,Vz B,Vz B sigma
c
c...third block
c...Bx A (rtn),Bx A sigma (rtn),Bx B (rtn),Bx B sigma (rtn)
c...By A (rtn),By A sigma (rtn),By B (rtn),By B sigma (rtn)
c...Bz A (rtn),Bz A sigma (rtn),Bz B (rtn),Bz B sigma (rtn)
c...Vx A (rtn),Vx A sigma (rtn),Vx B (rtn),Vx B sigma (rtn)
c...Vy A (rtn),Vy A sigma (rtn),Vy B (rtn),Vy B sigma (rtn)
c...Vz A (rtn),Vz A sigma (rtn),Vz B (rtn),Vz B sigma (rtn)
c.. Btrace A (rtn),Btrace B (rtn),Btrace A sigma (rtn),Btrace B sigma (rtn)
c...Vtrace A (rtn),Vtrace B (rtn),Vtrace As igma (rtn),Vtrace B sigma (rtn)
c...Vx mean, Vx Variance, Bx mean, Bx Variance
c...Vy mean, Vy variance, By mean, By variance
c...Vz mean, Vz variance, Bz mean, Bz variance
c...Vx variance**2, Vy variance**2, Vz_variance**2
c...Bx variance**2, By variance**2, Bz_variance**2
c
	write(69,*) Bx_mean,By_mean,Bz_mean
	write(18,900) iy1, t1, 
     .	       mean(4),var2(4),mean(1),var2(1),
     .		mean(2),var2(2),mean(3),var2(3), 
     .		mean(15),var2(15),  
     .		mean(11),var2(11),mean(10),var2(10), 
     .		mean(13),var2(13),mean(14),var2(14), 
     .		mean(8),var2(8),mean(9),var2(9),
     .	       mean(5),var2(5),mean(6),var2(6), 
     .	       mean(7),var2(7),Bvar,sigma_Bvar, 
     .	       var2(8)**2,var4(8),var2(5)**2, var4(5), 
     .         var2(6)**2,var4(6),var2(7)**2, var4(7), 
     .	       var2(4)**2,var4(4),var2(1)**2, var4(1), 
     .         var2(2)**2,var4(2),var2(3)**2, var4(3)

c
        write(18,901) Vvar,sigma_Vvar,  
     .		  VdotB_aver,VdotB_aver_var,VdotB_inst,VdotB_inst_var,  
     .		  Baniso,Vaniso,
     .            Btrace_A,Btrace_B,Btrace_sigmaA,Btrace_sigmaB,
     .            Vtrace_A,Vtrace_B,Vtrace_sigmaA,Vtrace_sigmaB,
     .	          Tp_A,Tp_sigmaA,Tp_B,Tp_sigmaB,
     .            Np_A,Np_sigmaA,Np_B,Np_sigmaB,
     .            Bx_A,Bx_sigmaA,Bx_B,Bx_sigmaB,
     .		  By_A,By_sigmaA,By_B,By_sigmaB,
     .            Bz_A,Bz_sigmaA,Bz_B,Bz_sigmaB,
     .            Vx_A,Vx_sigmaA,Vx_B,Vx_sigmaB,
     .	          Vy_A,Vy_sigmaA,Vy_B,Vy_sigmaB,
     .            Vz_A,Vz_sigmaA,Vz_B,Vz_sigmaB
     
        write(18,902)         Bx_A_rtn,Bx_sigmaA_rtn,Bx_B_rtn,Bx_sigmaB_rtn,
     .            By_A_rtn,By_sigmaA_rtn,By_B_rtn,By_sigmaB_rtn,
     .            Bz_A_rtn,Bz_sigmaA_rtn,Bz_B_rtn,Bz_sigmaB_rtn,
     .            Vx_A_rtn,Vx_sigmaA_rtn,Vx_B_rtn,Vx_sigmaB_rtn,
     .            Vy_A_rtn,Vy_sigmaA_rtn,Vy_B_rtn,Vy_sigmaB_rtn,
     .	          Vz_A_rtn,Vz_sigmaA_rtn,Vz_B_rtn,Vz_sigmaB_rtn,
     .	          Btrace_A_rtn,Btrace_B_rtn,Btrace_sigmaA_rtn,Btrace_sigmaB_rtn,
     .     	  Vtrace_A_rtn,Vtrace_B_rtn,Vtrace_sigmaA_rtn,Vtrace_sigmaB_rtn,
     .            Vx_mean,Vx_mean2,Bx_mean,Bx_mean2,
     .            Vy_mean,Vy_mean2,By_mean,By_mean2,
     .            Vz_mean,Vz_mean2,Bz_mean,Bz_mean2,
     .            Vx_mean4,Vy_mean4,Vz_mean4,
     .            Bx_mean4,By_mean4,Bz_mean4
c
 900	format(1x,i4,1x,f7.3,                                                     
     .	    /3x,2(1x,1pe12.4,1x,1pe12.4),
     .	    /3x,2(1x,1pe12.4,1x,1pe12.4),
     .	    /3x,1(1x,1pe12.4,1x,1pe12.4),	
     . 	    /3x,2(1x,1pe12.4,1x,1pe12.4),
     .      /3x,2(1x,1pe12.4,1x,1pe12.4), 
     .	    /3x,2(1x,1pe12.4,1x,1pe12.4), 
     .      /3x,2(1x,1pe12.4,1x,1pe12.4),
     .      /3x,2(1x,1pe12.4,1x,1pe12.4),  
     .	    /3x,2(1x,1pe12.4,1x,1pe12.4),	
     .      /3x,2(1x,1pe12.4,1x,1pe12.4),
     .      /3x,2(1x,1pe12.4,1x,1pe12.4), 
     .      /3x,2(1x,1pe12.4,1x,1pe12.4))
 

 901	format( 3x,1(1x,1pe12.4,1x,1pe12.4), 
     .      /3x,2(1x,1pe12.4,1x,1pe12.4), 
     .      /3x,1(1x,1pe12.4,1x,1pe12.4), 
     .      /3x,2(1x,1pe12.4,1x,1pe12.4),       
     .      /3x,2(1x,1pe12.4,1x,1pe12.4), 
     .      /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .	    /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4))
     
 902	format(3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .	    /3x,4(1x,1pe12.4),
     .      /3x,4(1x,1pe12.4),
     .      /3x,2(1x,1pe12.4,1x,1pe12.4), 
     .      /3x,2(1x,1pe12.4,1x,1pe12.4), 
     .      /3x,2(1x,1pe12.4,1x,1pe12.4), 
     .      /3x,3(1x,1pe12.4), 
     .      /3x,3(1x,1pe12.4))
c

	write(56,903) t2,Btrace_B
	write(57,903) t2,Btrace_sigmaB
	write(58,903) t2,Btrace_A
	write(59,903) t2,Btrace_sigmaA
 903	format(e14.5,3x,e14.5)
c.. get the next hour:
c	write(56,*) iy2,t2
	iy1 = iy2
	t1 = t2
	t2 = t2 + one_hour
	iy2 = iy2
	if(((iy1/4)*4 .ne. iy1) .and. (t2 .gt. 366.)) then
	    iy2 = iy2 + 1
	    t2 = t2 - 366.
	endif
	if(((iy1/4)*4 .eq. iy1) .and. (t2 .gt. 367.)) then
	    iy2 = iy2 + 1
	    t2 = t2 - 367.
	endif
	if(iy1 .gt. iyear_stop) go to 1000
	if((iy1 .eq. iyear_stop) .and. (t1 .gt. time_stop)) go to 1000
	go to 200
c
 1000	continue
	close(unit=14)
	close(unit=15)
	close(unit=18)
c	close(unit=90)
c	close(unit=91)
c	close(unit=92)
c	close(unit=93)
c	close(unit=94) 
c	close(unit=95)
c	close(unit=96)
c	close(unit=97)
c        close(unit=98)
c        close(unit=99)
c
	stop
	end
c
c
c
	subroutine gethour(iflag,iwhy,npts,Brms)
c
c.. data declarations for data input:
	integer*4 npts 
	integer*4 iyear, icnt_mag(60)
	real*8 time, Np(60), Tp(60), alpha(60), Vp(60), Vr(60), Vt(60), Vn(60)
	real*8 Br(60), Bt(60), Bn(60), Bmag(60), zlambda(60), delta_dat(60), Brms(60)
	real*8 Xsc(60), Ysc(60), Zsc(60)
	real*8 times
	common /hourly/ time, Np, Tp, alpha, Vp, Vr, Vt, Vn,
     .		Br, Bt, Bn, Bmag, zlambda, delta_dat, Xsc, Ysc, Zsc, 
     .		icnt_mag,times,iyear
c
c.. data declarations for data selection:
	integer*4 iyear_start, iyear_stop
	real*8 time_start, time_stop
	character*40 merged_files_listing
	common /timeselect/ iyear_start, iyear_stop, time_start, time_stop, 
     .		merged_files_listing
c
c.. data declarations for housekeeping variables:
	integer*4 iy1, iy2
	real*8 bad, pi, t1, t2
	common /house/ bad, pi, t1, t2, iy1, iy2
c
c.. passed variables:
	integer*4 iflag, iwhy
c
	iflag = 0
	iwhy = 0
c
	ii = 1
  200 continue
	  read(15,*,end=300) iyear, time, 
     .			Np(ii), Tp(ii), alpha(ii), Vp(ii), Vr(ii), Vt(ii), Vn(ii), 
     .			Br(ii), Bt(ii), Bn(ii), Bmag(ii), zlambda(ii), delta_dat(ii), Brms(ii), 
     .			Xsc(ii), Ysc(ii), Zsc(ii), icnt_mag(ii)
c
c.. exclude time out-of-bounds:
	  if(iyear .lt. iy1) go to 200
	  if((iyear .eq. iy1) .and. (time .lt. t1)) go to 200
	  if(iyear .gt. iy2) go to 500
	  if((iyear .eq. iy2) .and. (time .gt. t2)) go to 500
	  ii = ii + 1
	  go to 200
c
  300	  continue			! span datasets and flag out incomplete hour of data
	  call getfile(iflag,iwhy)
	  if(iflag .eq. 1) then
		go to 200
	    else
		return
	  endif
c..end of implied do loop for hourly construct
c
  500	continue			! return hour of data for analysis
	npts = ii - 1
	iflag = 1
	return
	end
c
c
c
	subroutine getfile(iflag,iwhy)
c
c.. data declarations for data selection:
	integer*4 iyear_start, iyear_stop
	real*8 time_start, time_stop
	character*40 merged_files_listing
	common /timeselect/ iyear_start, iyear_stop, time_start, time_stop, 
     .		merged_files_listing
c
c.. data file variables:
	character*3 aday1, aday2
	character*4 ayear1, ayear2
	character*60 aline
	character*60 aswepam_filename
	character*80 idget
	character*110 swepam_file
c
	iflag = 0					! reset to +/- 1
	iwhy = 0
c
c.. get MAG and SWEPAM data from merged files:
 100	continue
	read(14,110,end=300) aswepam_filename
 110	format(a60)
	ayear1 = aswepam_filename(21:24)
	read(ayear1,120) iyear1
 120	format(i4)
	ayear2 = aswepam_filename(30:33)
	read(ayear2,120) iyear2
	aday1 = aswepam_filename(26:28)
	read(aday1,130) iday1
 130	format(i3)
	aday2 = aswepam_filename(35:37)
	read(aday2,130) iday2
c.. select time intervals for processing:
	if(iyear2 .lt. iyear_start) go to 100
	if(iyear1 .gt. iyear_stop) go to 400
	if((iyear2 .eq. iyear_start) .and. (iday2 .lt. time_start)) go to 100
	if((iyear1 .eq. iyear_stop) .and. (iday1 .gt. time_stop)) go to 400
c	
	swepam_file = '/home/students/tessein/merged_ace_data/'//aswepam_filename
c	swepam_file = 'briaxa$dkb1200:[users.chuck.acelevel2data.swepam]'//aswepam_filename
	write(89,*) swepam_file,aswepam_filename
	open(unit=15, file=swepam_file, form='formatted', status='old')
	do 215 iskip=1,22
	  read(15,205) aline
  205	  format(a60)
	  write(6,210) aline
  210	  format(' ',a60)
  215	continue
c.. file is opened and header read - get back to work:
	iflag = 1
	return
c
 300	continue
	write(6,310)
 310	format('   End of processing interval - stop!')
	iflag = -1
	iwhy = 1
	return
 400	continue
	write(6,410)
 410	format('   Ran out of data files - stop!')
	iflag = -1
	iwhy = -1
	return
	end
c
c
      subroutine meanvar(y,nptss, ymean, yvar2, yvar4)
c...Computes means and variances    
      integer*4 nptss
      real*8 y(60), ymean,  yvar2,yvar4
      real*8 bad
c
      bad = -9999.9
      ymean = 0.
      yvar2 = 0.
      yvar4 = 0.
      icnt = 0
      do 100 i=1,nptss
	  if(abs(y(i)-bad) .gt. 1.) then
		ymean = ymean + y(i)
		yvar2 = yvar2 + y(i)**2
		yvar4 = yvar4 + y(i)**4
		icnt = icnt + 1
	  endif
  100	continue
	if(icnt .gt. 1) then
	    zcnt = float(icnt)
	    ymean = ymean / zcnt
	    yvar2 = yvar2 / zcnt
	    yvar4 = yvar4 / zcnt
	    yvar2 = yvar2 - ymean**2
	    yvar4 = yvar4 - yvar2**2
	    yvar2 = sqrt(yvar2)
	    yvar4 = sqrt(yvar4)
	  else
	    ymean = bad
	    yvar2 = bad
	    yvar4 = bad
	endif
c

	end

c
c
c---Routine MNFLD----------------------------------------------------
c
c  Activate using:
c	do 100 i=1,npts
c	  Bx(i) = Br(i)
c 100	continue
c	call mnfld(Bx,...)
c
c--------------------------------------------------------------------
c
       subroutine mnfld(npts,Baniso,Vaniso)
c........
c  data declaration for data input:
      integer*4 nser, npts, nhist, jtype, NDAT_b,ndat_v,ndat
        real*8 Bx(60), By(60), Bz(60)
        real*8 Vx(60), Vy(60), Vz(60)
	real*8 Baniso,Vaniso
	real*8 rm(3,3)
	real*8 bad, pi
        real*8 xsum, ysum, zsum, Vxm, Vym, Vzm, Bxm, Bym, Bzm 
	real*8 BxAv, ByAv, BzAv, BrSqr, BtSqr, Bnsqr, Bav, Blo, Bxy, Byz, Bxz, yr, ymag
	real*8 denom, costheta
	real*8 Bz_mean, Bz_mean2, Bz_mean4, Bx_mean, Bx_mean2, Bx_mean4, By_mean, By_mean2, By_mean4
	real*8 Vz_mean, Vz_mean2, Vz_mean4, Vx_mean, Vx_mean2, Vx_mean4, Vy_mean, Vy_mean2, Vy_mean4 
	integer*4 icnt_Bx, icnt_By, icnt_Bz, icnt_Vx, icnt_Vy, icnt_Vz
	integer*4 iyear_start,iyear_stop,i,j
	common /SWEPAM/ Bx, By, Bz, Vx, Vy, Vz
c
	common /mnfld_data/ Bx_mean,Bx_mean2,Bx_mean4,By_mean,By_mean2,By_mean4,Bz_mean,Bz_mean2,Bz_mean4,
     .			    Vx_mean,Vx_mean2,Vx_mean4,Vy_mean,Vy_mean2,Vy_mean4,Vz_mean,Vz_mean2,Vz_mean4
c
	real*8 time_start,time_stop
	character*40 merged_files_listing 
	common /timeselect/ iyear_start, iyear_stop, time_start, time_stop, 
     .	      	merged_files_listing
c
      bad = -9999.9
      pi = acos(-1.0)
      ndat_b = 0
      ndat_v = 0
      ndat = 0
      xsum = 0.0
      ysum = 0.0
      zsum = 0.0
c
c... compute mean vector:
        do 4 i=1,npts
	    if((abs(Bx(i) - bad) .gt. 0.1) .and. 
     .       (abs(By(i) - bad) .gt. 0.1) .and. 
     .       (abs(Bz(i) - bad) .gt. 0.1)) then
    		xsum = xsum + Bx(i)
		ysum = ysum + By(i)
		zsum = zsum + Bz(i)
		ndat = ndat + 1
	    endif
 4      continue
        if(ndat .gt. 1) then
	      BxAv = xsum/ndat
            BrSqr = BxAv**2
	      ByAv = ysum/ndat
            BtSqr = ByAv**2
	      BzAv = zsum/ndat
            BnSqr = BzAv**2
	      Bav = sqrt( BxAv**2 + ByAv**2 + BzAv**2 )
	      Blo = 0.01*Bav
	      Bxy = sqrt(BxAv**2+ByAv**2)
	      Byz = sqrt(ByAv**2+BzAv**2)
	      Bxz = SQRT(BxAv**2+BzAv**2)
            yr = BtSqr + BnSqr
            ymag = sqrt(BtSqr**2 + BtSqr*BrSqr + 2*BtSqr*BnSqr + BnSqr*BrSqr + BnSqr**2)      
	      costheta = yr/ymag
          else
	      BxAv = bad
	      BrSqr = bad
	      ByAv = bad
 	      BtSqr = bad
	      BzAv = bad
	      BnSqr = bad
	      Bav = bad
	      Blo = bad
	      Bxy = bad
	      Byz = bad
	      yr = bad
	      ymag = bad
	      costheta = bad
	  endif
c
c Added code for IPLANE=4 (Bieber et al. rotation)
c RJL, 19970518
c
	if(ndat .gt. 1) then
	    denom = sqrt(Bav**2 * (Bav**2 - BxAv**2))
	    rm(1,1) = 0.0
	    rm(1,2) = BzAv/Byz
	    rm(1,3) = -ByAv/Byz
	    rm(2,1) = -(Bav**2 - BxAv**2)/denom
	    rm(2,2) = BxAv*ByAv /denom
	    rm(2,3) = BxAv*BzAv /denom
	    rm(3,1) = BxAv/Bav
	    rm(3,2) = ByAv/Bav
	    rm(3,3) = BzAv/Bav
	else
	    Do 6  i=1,npts
	       Vx(i) = bad
	       Vy(i) = bad
	       Vz(i) = bad
 6	    continue
        endif
c
  	
c	
c	do 10 i=1,npts
c	  if((abs(Bx(i) - bad) .gt. 0.1) .and.
c     .     (abs(By(i) - bad) .gt. 0.1) .and. 
c     .     (abs(Bz(i) - bad) .gt. 0.1)) then
c	       bx(i) = bx(i) - bxav
c	       by(i) = by(i) - byav
c	       bz(i) = bz(i) - bzav
c        endif
c   10	continue
c

	vxav = 0.
	vyav = 0.
	vzav = 0.
	do 15 i=1,npts
	  if((abs(Vx(i) - bad) .gt. 0.1) .and.
     .   (abs(Vy(i) - bad) .gt. 0.1) .and.
     .   (abs(Vz(i) - bad) .gt. 0.1)) then
    		vxav = vxav + Vx(i)
		vyav = vyav + Vy(i)
		vzav = vzav + Vz(i)
		ndat = ndat + 1 
	  endif
   15	continue
	if(ndat .eq. 0) then
	   VxAv = bad
	   VyAv = bad
	   VzAv = bad
	  else 
	   VxAv = vxav / float(ndat)
	   VyAv = vyav / float(ndat)
	   VzAv = vzav / float(ndat)
	endif
c	do 20 i=1,npts
c	  if((abs(Vx(i) - bad) .gt. 0.1) .and. 
c     .     (abs(Vy(i) - bad) .gt. 0.1) .and.
c     .     (abs(Vz(i) - bad) .gt. 0.1)) then
c    		Vx(i) = Vx(i) - vxav
c		Vy(i) = Vy(i) - vyav
c		Vz(i) = Vz(i) - vzav
c	  endif
c   20	continue
c
c Rotate vector
c
	DO 30 I=1,npts
	  if ((abs(Bx(i)-bad) .lt. 0.1) .or. 
     .      (abs(Bx(i)-bad) .lt. 0.1) .or.
     .	(abs(Bx(i)-bad) .lt. 0.1) .or. 
     .	(ndat .le. 1)) then
		Bxm = bad
		Bym = bad
		Bzm = bad
	    else
		Bxm = rm(1,1)*Bx(i)+rm(1,2)*By(i)+rm(1,3)*Bz(i)
		Bym = rm(2,1)*Bx(i)+rm(2,2)*By(i)+rm(2,3)*Bz(i)
		Bzm = rm(3,1)*Bx(i)+rm(3,2)*By(i)+rm(3,3)*Bz(i)
        endif
	  Bx(i) = Bxm
	  By(i) = Bym
	  Bz(i) = Bzm
   30	CONTINUE
c
	DO 40 I=1,NPTS
	  if ((abs(Vx(i)-bad) .lt. 0.1) .or. 
     .      (abs(Vy(i)-bad) .lt. 0.1) .or.
     .	(abs(Vz(i)-bad) .lt. 0.1) .or. 
     .	(ndat .le. 1)) then
		Vxm = bad
		Vym = bad
		Vzm = bad 
	    else
		Vxm = rm(1,1)*Vx(i)+rm(1,2)*Vy(i)+rm(1,3)*Vz(i)
		Vym = rm(2,1)*Vx(i)+rm(2,2)*Vy(i)+rm(2,3)*Vz(i)
		Vzm = rm(3,1)*Vx(i)+rm(3,2)*Vy(i)+rm(3,3)*Vz(i)
        endif
	  Vx(i) = Vxm
	  Vy(i) = Vym
	  Vz(i) = Vzm
 40   CONTINUE
c.....Initialize statistics variables to 0
	Bx_mean = 0.
	Bx_mean2 = 0.
	Bx_mean4 = 0.
	By_mean = 0.
	By_mean2 = 0.
	By_mean4 = 0.
	Bz_mean = 0.
	Bz_mean2 = 0.
	Bz_mean4 = 0.
	Vx_mean = 0.
	Vx_mean2 = 0.
	Vx_mean4 = 0.
	Vy_mean = 0.
	Vy_mean2 = 0.
	Vy_mean4 = 0.
	Vz_mean = 0.
	Vz_mean2 = 0.
	Vz_mean4 = 0.
c... Initialize counting variables to 0
	icnt_Bx = 0
	icnt_By = 0
	icnt_Bz = 0
	icnt_Vx = 0
	icnt_Vy = 0
	icnt_Vz = 0
c.....Running 1-hour sums
 	do 90 i=1,npts
	  if(abs(Bx(i)-bad) .gt. 1.) then
		Bx_mean = Bx_mean + Bx(i)
		Bx_mean2 = Bx_mean2 + Bx(i)**2
		Bx_mean4 = Bx_mean4 + Bx(i)**4
		icnt_Bx = icnt_Bx + 1
	  endif
	  if(abs(By(i)-bad) .gt. 1.) then
	      By_mean = By_mean + By(i)
	      By_mean2 = By_mean2 + By(i)**2
	      By_mean4 = By_mean4 + By(i)**4
	      icnt_By = icnt_By + 1
	  endif
	  if(abs(Bz(i)-bad) .gt. 1.) then
	      Bz_mean = Bz_mean + Bz(i)
	      Bz_mean2 = Bz_mean2 + Bz(i)**2
	      Bz_mean4 = Bz_mean4 + Bz(i)**4
	      icnt_Bz = icnt_Bz + 1
	  endif
  	  if(abs(Vx(i)-bad) .gt. 1.) then
  		Vx_mean = Vx_mean + Vx(i)
  		Vx_mean2 = Vx_mean2 + Vx(i)**2
  		Vx_mean4 = Vx_mean4 + Vx(i)**4
  		icnt_Vx = icnt_Vx + 1
  	  endif
  	  if(abs(Vy(i)-bad) .gt. 1.) then
  		Vy_mean = Vy_mean + Vy(i)
  		Vy_mean2 = Vy_mean2 + Vy(i)**2
  		Vy_mean4 = Vy_mean4 + Vy(i)**4
  		icnt_Vy = icnt_Vy + 1
  	  endif
  	  if(abs(Vz(i)-bad) .gt. 1.) then
  		Vz_mean = Vz_mean + Vz(i)
  		Vz_mean2 = Vz_mean2 + Vz(i)**2
  		Vz_mean4 = Vz_mean4 + Vz(i)**4
		icnt_Vz = icnt_Vz + 1
	  endif
 90	continue
c.......Compute means and variances
	  if(icnt_Bx .gt. 1) then
            Bx_mean = Bx_mean / float(icnt_Bx)
		Bx_mean2 = Bx_mean2 / float(icnt_Bx)
		Bx_mean4 = Bx_mean4 / float(icnt_Bx)
		Bx_mean2 = Bx_mean2 - Bx_mean**2
		Bx_mean4 = Bx_mean4 - Bx_mean2**2
		Bx_mean2 = sqrt(Bx_mean2)
		Bx_mean4 = sqrt(Bx_mean4)
	  else
		Bx_mean = bad
		Bx_mean2 = bad
		Bx_mean4 = bad
	  endif
	  if(icnt_By .gt. 1) then	
		By_mean = By_mean / float(icnt_By)
		By_mean2 = By_mean2 / float(icnt_By)
		By_mean4 = By_mean4 / float(icnt_By)
		By_mean2 = By_mean2 - By_mean**2
		By_mean4 = By_mean4 - By_mean2**2
		By_mean2 = sqrt(By_mean2)
		By_mean4 = sqrt(By_mean4)	
	  else
		By_mean = bad
		By_mean2 = bad
		By_mean4 = bad
	  endif
	  if(icnt_Bz .gt. 1) then	
		Bz_mean = Bz_mean / float(icnt_Bz)
	Bz_mean2 = Bz_mean2 / float(icnt_Bz)
		Bz_mean4 = Bz_mean4 / float(icnt_Bz)
		Bz_mean2 = Bz_mean2 - Bz_mean**2
		Bz_mean4 = Bz_mean4 - Bz_mean2**2
		Bz_mean2 = sqrt(Bz_mean2)
		Bz_mean4 = sqrt(Bz_mean4)
	  else
	      Bz_mean = bad
	      Bz_mean2 = bad
	      Bz_mean4 = bad
	  endif
	  if(icnt_Vx .gt. 1) then
	      Vx_mean = Vx_mean / float(icnt_Vx)
	      Vx_mean2 = Vx_mean2 / float(icnt_Vx)
	      Vx_mean4 = Vx_mean4 / float(icnt_Vx)
	      Vx_mean2 = Vx_mean2 - Vx_mean**2
	      Vx_mean4 = Vx_mean4 - Vx_mean2**2
	      Vx_mean2 = sqrt(Vx_mean2)
	      Vx_mean4 = sqrt(abs(Vx_mean4))	
	  else
	      Vx_mean = bad
	      Vx_mean2 = bad
	      Vx_mean4 = bad
	  endif
	  if(icnt_Vy .gt. 2) then
	      Vy_mean = Vy_mean / float(icnt_Vy)
	      Vy_mean2 = Vy_mean2 / float(icnt_Vy)
 	      Vy_mean4 = Vy_mean4 / float(icnt_Vy)
            Vy_mean2 = Vy_mean2 - Vy_mean**2
            Vy_mean4 = Vy_mean4 - Vy_mean2**2
            Vy_mean2 = sqrt(Vy_mean2)
	      Vy_mean4 = sqrt(Vy_mean4)	
	  else
	      Vy_mean = bad
	      Vy_mean2 = bad
	      Vy_mean4 = bad
	  endif
	  if(icnt_Vz .gt. 2) then
	      Vz_mean = Vz_mean / float(icnt_Vz)
	      Vz_mean2 = Vz_mean2 / float(icnt_Vz)
	      Vz_mean4 = Vz_mean4 / float(icnt_Vz)
	      Vz_mean2 = Vz_mean2 - Vz_mean**2
	      Vz_mean4 = Vz_mean4 - Vz_mean2**2
	      Vz_mean2 = sqrt(Vz_mean2)
	      Vz_mean4 = sqrt(Vz_mean4)
	   else
	      Vz_mean = bad
	      Vz_mean2 = bad 
	      Vz_mean4 = bad
         endif
c.....Compute B and V anisotropy
	if(((Bx_mean2-bad) .gt. 1.) .and. ((By_mean2-bad) .gt. 1.) .and. ((Bz_mean2-bad) .gt. 1.)) then
	    Baniso = (Bx_mean2 + By_mean2) / Bz_mean2
	  else
	    Baniso = bad
	endif
	if(((Vx_mean2-bad) .gt. 1.) .and. ((Vy_mean2-bad) .gt. 1.) .and. 
     .     ((Vz_mean2-bad) .gt. 1.) .and. (Vz_mean2 .gt. 0.)) then    
		Vaniso = (Vx_mean2 + Vy_mean2) / Vz_mean2
	else
	    	Vaniso = bad
	 endif
c
	return
      end 
c
c
c....structure functions routine

	subroutine structurefunc(npts)
c....Subroutine to compute structure functions
c
c...J.A.T. -- 7/18/06
c
c...Variable Declaration
	real*8 Bx(60), By(60), Bz(60), Vx(60), Vy(60), Vz(60), Np(60), Tp(60)
	real*8 S2_Bx(60), S2_By(60), S2_Bz(60), S2_Vx(60), S2_Vy(60), S2_Vz(60), S2_Np(60), S2_Tp(60)
	real*8 S2_Bx2(60), S2_By2(60), S2_Bz2(60), S2_Vx2(60), S2_Vy2(60), S2_Vz2(60), S2_Np2(60), S2_Tp2(60)
 	real*8 Bx_StrucTemp, By_StrucTemp, Bz_StrucTemp, Vx_StrucTemp, Vy_StrucTemp, Vz_StrucTemp
  	real*8 Np_StrucTemp, Tp_StrucTemp, Btrace(60), Vtrace(60),Btrace2(60),Vtrace2(60)
	real*8 Bx_StrucTemp2,By_StrucTemp2,Bz_StrucTemp2,Vx_StrucTemp2,Vy_StrucTemp2,Vz_StrucTemp2
	real*8 Np_StrucTemp2,Tp_StrucTemp2
	real*8 Btrace_StrucTemp,Btrace_StrucTemp2,Vtrace_StrucTemp,Vtrace_StrucTemp2
	real*8 bad,temp,temp2,delta_dat(60)
	real*8 Tp_A,Tp_sigmaA,Tp_B,Tp_sigmaB,Np_A,Np_sigmaA,Np_B,Np_sigmaB
        real*8 Bx_A,Bx_sigmaA,Bx_B,Bx_sigmaB,By_A,By_sigmaA,By_B,By_sigmaB
        real*8 Bz_A,Bz_sigmaA,Bz_B,Bz_sigmaB,Vx_A,Vx_sigmaA,Vx_B,Vx_sigmaB
     	real*8 Vy_A,Vy_sigmaA,Vy_B,Vy_sigmaB,Vz_A,Vz_sigmaA,Vz_B,Vz_sigmaB
     	real*8 Btrace_A, Btrace_B,Btrace_sigmaA,Btrace_sigmaB
     	real*8 Vtrace_A, Vtrace_B,Vtrace_sigmaA,Vtrace_sigmaB
     	real*8 sigma_BxStruc(60),sigma_ByStruc(60),sigma_BzStruc(60),sigma_BtraceStruc(60),sigma_VtraceStruc(60)
     	real*8 sigma_VxStruc(60),sigma_VyStruc(60),sigma_VzStruc(60),sigma_NpStruc(60),sigma_TpStruc(60)
	integer*4 icnt_VxStruc, icnt_VyStruc, icnt_VzStruc, icnt_NpStruc, icnt_TpStruc, icnt_Bstruc 
	integer*4 N, k, i, j,iyear,icnt_mag(60),icnt_VtraceStruc
c.. native data declarations/common blocks:      
	common /SWEPAM/ Bx, By, Bz, Vx, Vy, Vz
	common /structure/ S2_Tp, S2_Np, S2_Bx, S2_By, S2_Bz, S2_Vx, S2_Vy, S2_Vz,k,N,Btrace,Vtrace,
     .                     S2_Tp2, S2_Np2, S2_Bx2, S2_By2, S2_Bz2, S2_Vx2, S2_Vy2, S2_Vz2,sigma_VxStruc,
     .                     sigma_VyStruc,sigma_VzStruc,sigma_NpStruc,sigma_TpStruc,sigma_VtraceStruc,
     .                     sigma_BxStruc,sigma_ByStruc,sigma_BzStruc,sigma_BtraceStruc
	real*8 time,alpha(60),Vp(60),Vr(60),Vt(60),Vn(60)
	real*8 zlambda(60),delta(60,15),Br(60),Bt(60),Bn(60),Bmag(60)
	real*8 Xsc(60),Ysc(60),Zsc(60),times
c
	common /hourly/ time, Np, Tp, alpha, Vp, Vr, Vt, Vn,
     .		Br, Bt, Bn, Bmag, zlambda, delta_dat, Xsc, Ysc, Zsc, 
     .		icnt_mag,times,iyear
c
	real*8 delta_fit(10),A_fit,B_fit,variance,sigmaA,zn,sigma_zn,Zamp,sigma_Zamp,sigmaB
c
	common /linefit_var/ delta_fit,A_fit,B_fit,variance,sigmaA,sigmaB,zn,sigma_zn,Zamp,sigma_Zamp
	common /passfit/ Tp_A,Tp_sigmaA,Tp_B,Tp_sigmaB,Np_A,Np_sigmaA,Np_B,Np_sigmaB,
     .                   Bx_A,Bx_sigmaA,Bx_B,Bx_sigmaB,By_A,By_sigmaA,By_B,By_sigmaB,
     .                   Bz_A,Bz_sigmaA,Bz_B,Bz_sigmaB,Vx_A,Vx_sigmaA,Vx_B,Vx_sigmaB,
     .                   Vy_A,Vy_sigmaA,Vy_B,Vy_sigmaB,Vz_A,Vz_sigmaA,Vz_B,Vz_sigmaB,
     .                   Btrace_A,Btrace_sigmaA,Btrace_B,Btrace_sigmaB,
     .                   Vtrace_A,Vtrace_sigmaA,Vtrace_B,Vtrace_sigmaB
c
	
c....initialize variables
	N = 57 
	bad = -9999.9
	do 10 i=1,60
	  S2_Bx(i) = bad
	  S2_By(i) = bad
	  S2_Bz(i) = bad
	  S2_Vx(i) = bad
	  S2_Vy(i) = bad
	  S2_Vz(i) = bad
	  S2_Np(i) = bad
	  S2_Tp(i) = bad
	  S2_Bx2(i) = bad
	  S2_By2(i) = bad
	  S2_Bz2(i) = bad
	  S2_Vx2(i) = bad
	  S2_Vy2(i) = bad
	  S2_Vz2(i) = bad
	  S2_Np2(i) = bad
	  S2_Tp2(i) = bad
   10	continue
c                 
c....Begin computing structure functions, k is the lag
	Do 500 k = 1,N/5
 	 icnt_Bstruc  = 0
	 icnt_VxStruc = 0
	 icnt_VyStruc = 0
	 icnt_VzStruc = 0
	 icnt_NpStruc = 0
	 icnt_TpStruc = 0
	 icnt_VtraceStruc = 0
	 i = 0
	 j = 0
	 Bx_StrucTemp = 0.
	 By_StrucTemp = 0. 
	 Bz_StrucTemp = 0. 
	 Vx_StrucTemp = 0. 
	 Vy_StrucTemp = 0. 
	 Vz_StrucTemp = 0. 
	 Np_StrucTemp = 0. 
	 Tp_StrucTemp = 0.
	 Bx_StrucTemp2 = 0.
	 By_StrucTemp2 = 0. 
	 Bz_StrucTemp2 = 0. 
	 Vx_StrucTemp2 = 0. 
	 Vy_StrucTemp2 = 0. 
	 Vz_StrucTemp2 = 0. 
	 Np_StrucTemp2 = 0. 
	 Tp_StrucTemp2 = 0.
	 Btrace_StrucTemp = 0.
	 Btrace_StrucTemp2 = 0.
	 Vtrace_StrucTemp = 0.
	 Vtrace_StrucTemp2 = 0.
c	 
 100	 continue
 	 i = i + 1
 	 j = i + k
c...Begin computing structure functions here
       if((abs(Bx(i)-bad) .gt. 1.) .and. (abs(Bx(j)-bad) .gt. 1.)) then
 	    Bx_StrucTemp = Bx_StrucTemp + (Bx(i) - Bx(j))**2
 	    Bx_StrucTemp2 = Bx_StrucTemp2 + (Bx(i) - Bx(j))**4
 	    icnt_Bstruc = icnt_Bstruc + 1
 	    By_StrucTemp = By_StrucTemp + (By(i) - By(j))**2
 	    By_StrucTemp2 = By_StrucTemp2 + (By(i) - By(j))**4
 	    Bz_StrucTemp = Bz_StrucTemp + (Bz(i) - Bz(j))**2
 	    Bz_StrucTemp2 = Bz_StrucTemp2 + (Bz(i) - Bz(j))**4
 	    temp = (Bx(i) - Bx(j))**2 + (By(i) - By(j))**2 + (Bz(i) - Bz(j))**2
 	    temp2 = temp**2
          Btrace_StrucTemp = Btrace_StrucTemp + temp
          Btrace_StrucTemp2 = Btrace_StrucTemp2 + temp2
	 endif 
   	 if((abs(Vx(i)-bad) .gt. 1.) .and. (abs(Vx(j)-bad) .gt. 1.)) then
   	    Vx_StrucTemp = Vx_StrucTemp + (Vx(i) - Vx(j))**2
   	    Vx_StrucTemp2 = Vx_StrucTemp2 + (Vx(i) - Vx(j))**4
   	    icnt_VxStruc = icnt_VxStruc + 1
       endif
       if((abs(Vy(i)-bad) .gt. 1.) .and. (abs(Vy(j)-bad) .gt. 1.)) then
     	    Vy_StrucTemp = Vy_StrucTemp + (Vy(i) - Vy(j))**2
     	    Vy_StrucTemp2 = Vy_StrucTemp2 + (Vy(i) - Vy(j))**4
     	    temp = (Vx(i) - Vx(j))**2 + (Vy(i) - Vy(j))**2 + (Vz(i) - Vz(j))**2
            temp2 = temp**2
            icnt_VyStruc = icnt_VyStruc + 1
	 endif   	
     	 if((abs(Vz(i)-bad) .gt. 1.) .and. (abs(Vz(j)-bad) .gt. 1.)) then
     	    Vz_StrucTemp = Vz_StrucTemp + (Vz(i) - Vz(j))**2
     	    Vz_StrucTemp2 = Vz_StrucTemp2 + (Vz(i) - Vz(j))**4
     	    icnt_VzStruc = icnt_VzStruc + 1    
	 endif
	 if((abs(Vx(i)-bad) .gt. 1.) .and. (abs(Vx(j)-bad) .gt. 1.) .and.
     .      (abs(Vy(i)-bad) .gt. 1.) .and. (abs(Vy(j)-bad) .gt. 1.) .and.
     .      (abs(Vz(i)-bad) .gt. 1.) .and. (abs(Vz(j)-bad) .gt. 1.)) then   
               temp = (Vx(i) - Vx(j))**2 + (Vy(i) - Vy(j))**2 + (Vz(i) - Vz(j))**2
               temp2 = temp**2
	       Vtrace_StrucTemp = Vtrace_StrucTemp + temp
               Vtrace_StrucTemp2 = Vtrace_StrucTemp2 + temp2
	       icnt_VtraceStruc = icnt_VtraceStruc + 1
	 endif
     	 if((abs(Np(i)-bad) .gt. 1.) .and. (abs(Np(j)-bad) .gt. 1.)) then
   	    Np_StrucTemp = Np_StrucTemp + (Np(i) - Np(j))**2
   	    Np_StrucTemp2 = Np_StrucTemp2 + (Np(i) - Np(j))**4
   	    icnt_NpStruc = icnt_NpStruc + 1
   	 endif   	
     	 if((abs(Tp(i)-bad) .gt. 1.) .and. (abs(Tp(j)-bad) .gt. 1.)) then
     	    Tp_StrucTemp = Tp_StrucTemp + (Tp(i) - Tp(j))**2
     	    Tp_StrucTemp2 =  Tp_StrucTemp2 + (Tp(i) - Tp(j))**4
     	    icnt_TpStruc = icnt_TpStruc + 1
	 endif
	 if(j .lt. N) goto 100
c
c...Finish computing structure functions and uncertainty for each
	 if(icnt_Bstruc .gt. 1) then
	     S2_Bx(k) = Bx_StrucTemp / float(icnt_Bstruc)
	     S2_Bx2(k) = Bx_StrucTemp2 / float(icnt_BStruc)
	     S2_By(k) = By_StrucTemp / float(icnt_Bstruc)
	     S2_By2(k) = By_StrucTemp2 / float(icnt_BStruc)
	     S2_Bz(k) = Bz_StrucTemp / float(icnt_Bstruc)
	     S2_Bz2(k) = Bz_StrucTemp2 / float(icnt_BStruc)
	     Btrace(k) = Btrace_StrucTemp / float(icnt_BStruc)
           Btrace2(k) = Btrace_StrucTemp2 / float(icnt_BStruc)
	     sigma_BxStruc(k) = sqrt((S2_Bx2(k) - (S2_Bx(k))**2)/float(icnt_BStruc))
	     sigma_ByStruc(k) = sqrt((S2_By2(k) - (S2_By(k))**2)/float(icnt_BStruc))
	     sigma_BzStruc(k) = sqrt((S2_Bz2(k) - (S2_Bz(k))**2)/float(icnt_BStruc))     
	     sigma_BtraceStruc(k) = sqrt((Btrace2(k) - (Btrace(k))**2)/float(icnt_BStruc))
	 else
	     S2_Bx(k) = bad
	     S2_By(k) = bad
	     S2_Bz(k) = bad
	     Btrace(k) = bad
	     Btrace2(k) = bad
	     S2_Bx2(k) = bad
	     S2_By2(k) = bad
 	     S2_Bz2(k) = bad
	     sigma_BxStruc(k) = bad 
	     sigma_ByStruc(k) = bad
	     sigma_BzStruc(k) = bad
	     sigma_BtraceStruc(k) = bad
	 endif
	 if(icnt_VxStruc .gt. 1) then
	 	S2_Vx(k) = Vx_StrucTemp / float(icnt_VxStruc)
		S2_Vx2(k) = Vx_StrucTemp2 / float(icnt_VxStruc)
		sigma_VxStruc(k) = sqrt((S2_Vx2(k) - (S2_Vx(k))**2)/float(icnt_VxStruc))
	 else
	      S2_Vx(k) = bad
	      S2_Vx2(k) = bad
	        sigma_VxStruc(k) = bad
	 endif
	 if(icnt_VyStruc .gt. 1) then
	 	S2_Vy(k) = Vy_StrucTemp / float(icnt_VyStruc)
	 	S2_Vy2(k) = Vy_StrucTemp2 / float(icnt_VyStruc)
	 	sigma_VyStruc(k) = sqrt((S2_Vy2(k) - (S2_Vy(k))**2)/float(icnt_VyStruc))
	 else
	 	S2_Vy(k) = bad
	 	S2_Vy2(k) = bad
	 	sigma_VyStruc(k) = bad
	 endif
	 if(icnt_VzStruc .gt. 1) then
	 	S2_Vz(k) = Vz_StrucTemp / float(icnt_VzStruc)
	 	S2_Vz2(k) = Vz_StrucTemp2 / float(icnt_VzStruc)
	 	sigma_VzStruc(k) = sqrt((S2_Vz2(k) - (S2_Vz(k))**2)/float(icnt_VzStruc))
	 else
	 	S2_Vz(k) = bad
	 	S2_Vz2(k) = bad
	 	sigma_VzStruc(k) = bad
	 endif
	 if(icnt_NpStruc .gt.1) then
	 	S2_Np(k) = Np_StrucTemp / float(icnt_NpStruc)
	 	S2_Np2(k) = Np_StrucTemp2 / float(icnt_NpStruc)
	 	sigma_NpStruc(k) = sqrt((S2_Np2(k) - (S2_Np(k))**2)/float(icnt_NpStruc))
       else
       	S2_Np(k) = bad
       	S2_Np2(k) = bad
       	sigma_NpStruc(k) = bad
       endif
       if(icnt_TpStruc .gt. 1) then
       	S2_Tp(k) = Tp_StrucTemp / float(icnt_TpStruc)
       	S2_Tp2(k) = Tp_StrucTemp2 / float(icnt_TpStruc)
       	sigma_TpStruc(k) = sqrt((S2_Tp2(k) - (S2_Tp(k))**2)/float(icnt_TpStruc))
	 else
	 	S2_Tp(k) = bad
	 	S2_Tp2(k) = bad
	 	sigma_TpStruc(k) = bad
	 endif
c...Compute trace
	 if((icnt_VxStruc .gt. 1) .and. (icnt_VyStruc .gt. 1) .and. (icnt_VzStruc .gt. 1)) then  
     	     Vtrace(k) = Vtrace_StrucTemp / float(icnt_VtraceStruc)
     	     Vtrace2(k) = Vtrace_StrucTemp2 / float(icnt_VtraceStruc)
     	     sigma_VtraceStruc(k) = sqrt((Vtrace2(k) - (Vtrace(k))**2)/float(icnt_VxStruc))
       else
	     Vtrace(k) = bad
	     Vtrace2(k) = bad
	     sigma_VtraceStruc(k) = bad
	 endif
 500 	continue
      k = k-1    !we can use k in the linefit
	call linefit
        return
	end



	 subroutine linefit 
c
c....Fit a regression to the structure function data
c
c Modified to compute linear fits to log/log scaling of the
c   observed dependence of P_perp/P_par upon proton beta.
c
c Uses least-squares routine pulled from Bevington.
c All data equally weighted, although this puts greater emphasis
c on the higher frequencies.
c
c   y = a + bx
c
c.... J.A.T. -- 7/26/06
c
	integer*4 npts, isum,i,count,k,N
	real*8 llog(2000),y_aver,taulog(60) 
	real*8 tau(60), tau_av
	real*8 tau_y_aver
	real*8 y_sqrd,x,y,x2,y2,xy,znpts
	real*8 sigma_zn, sigmaA, sigmaB, Zamp, sigma_Zamp, variance, A_fit, B_fit, delta_fit(10)
	real*8 sigma_A,sigma_B,zn
	real*8 bad, Btrace(60), Vtrace(60)
	real*8 tau_temp
	real*8 S2_Bx(60), S2_By(60), S2_Bz(60)
	real*8 S2_Vx(60), S2_Vy(60), S2_Vz(60), S2_Np(60), S2_Tp(60)
	real*8 Tp_A,Tp_sigmaA,Tp_B,Tp_sigmaB,Np_A,Np_sigmaA,Np_B,Np_sigmaB
        real*8 Bx_A,Bx_sigmaA,Bx_B,Bx_sigmaB,By_A,By_sigmaA,By_B,By_sigmaB
        real*8 Bz_A,Bz_sigmaA,Bz_B,Bz_sigmaB,Vx_A,Vx_sigmaA,Vx_B,Vx_sigmaB
     	real*8 Vy_A,Vy_sigmaA,Vy_B,Vy_sigmaB,Vz_A,Vz_sigmaA,Vz_B,Vz_sigmaB
        real*8 Btrace_A,Btrace_sigmaA,Btrace_B,Btrace_sigmaB
        real*8 Vtrace_A,Vtrace_sigmaA,Vtrace_B,Vtrace_sigmaB
        real*8 linefit_data(60,10),linefit_sigma(60,10),Sx(10),S(10),Sy(10),Sxx(10),Sxy(10)
        real*8 S2_Tp2(60),S2_Np2(60),S2_Bx2(60),S2_By2(60),S2_Bz2(60),S2_Vx2(60),S2_Vy2(60),S2_Vz2(60)
	real*8 sigma_VxStruc(60),sigma_VyStruc(60),sigma_VzStruc(60),sigma_NpStruc(60),sigma_TpStruc(60)
	real*8 sigma_BxStruc(60),sigma_ByStruc(60),sigma_BzStruc(60),sigma_VtraceStruc(60),sigma_BtraceStruc(60)
c
	common /structure/ S2_Tp, S2_Np, S2_Bx, S2_By, S2_Bz, S2_Vx, S2_Vy, S2_Vz,k,N,Btrace,Vtrace,
     .                     S2_Tp2, S2_Np2, S2_Bx2, S2_By2, S2_Bz2, S2_Vx2, S2_Vy2, S2_Vz2,sigma_VxStruc,
     .                     sigma_VyStruc,sigma_VzStruc,sigma_NpStruc,sigma_TpStruc,sigma_VtraceStruc,
     .                     sigma_BxStruc,sigma_ByStruc,sigma_BzStruc,sigma_BtraceStruc
	common /linefit_var/ delta_fit,A_fit,B_fit,variance,sigmaA,sigmaB,zn,sigma_zn,Zamp,sigma_Zamp
	common /linefit_stats/ x,y,x2,y2,xy,znpts
	common /passfit/ Tp_A,Tp_sigmaA,Tp_B,Tp_sigmaB,Np_A,Np_sigmaA,Np_B,Np_sigmaB,
     .                   Bx_A,Bx_sigmaA,Bx_B,Bx_sigmaB,By_A,By_sigmaA,By_B,By_sigmaB,
     .                   Bz_A,Bz_sigmaA,Bz_B,Bz_sigmaB,Vx_A,Vx_sigmaA,Vx_B,Vx_sigmaB,
     .                   Vy_A,Vy_sigmaA,Vy_B,Vy_sigmaB,Vz_A,Vz_sigmaA,Vz_B,Vz_sigmaB,
     .                   Btrace_A,Btrace_sigmaA,Btrace_B,Btrace_sigmaB,
     .                   Vtrace_A,Vtrace_sigmaA,Vtrace_B,Vtrace_sigmaB
c...initialize variables
	bad = -9999.9
	tau_temp = 64.	
c...Bring data into a square array linefit_data(60,10)
c...Columns: Tp,Np,Bx,By,Bz,Vx,Vy,Vz,Btrace,Vtrace
c
	do 10 i=1,k
	  linefit_data(i,1) = S2_Tp(i)
	  linefit_data(i,2) = S2_Np(i)
	  linefit_data(i,3) = S2_Bx(i)
	  linefit_data(i,4) = S2_By(i)
	  linefit_data(i,5) = S2_Bz(i)
	  linefit_data(i,6) = S2_Vx(i)
	  linefit_data(i,7) = S2_Vy(i)
	  linefit_data(i,8) = S2_Vz(i)
	  linefit_data(i,9) = Btrace(i)
	  linefit_data(i,10) = Vtrace(i)  
 10 	continue
c...Initialize arrays
	Do 1000 count=1,10
	 delta_fit(count) = 0.
	 S(count) = 0.
	 Sx(count) = 0.
	 Sy(count) = 0.
	 Sxx(count) = 0.
	 Sxy(count) = 0.
c
c...compute uncertainty
	isum = 0
	Do 77 i=1,k
	  if(abs(linefit_data(i,count)-bad) .gt. 1. ) then
	    isum = isum + 1
	    tau(i) = tau_temp * float(i)
	    taulog(i) = log10(tau(i))
 	    llog(i) = log10(linefit_data(i,count))
            if(count .eq. 1) linefit_sigma(i,count)=sigma_TpStruc(i)/(S2_Tp(i)*log(10.))  
	    if(count .eq. 2) linefit_sigma(i,count)=sigma_NpStruc(i)/(S2_Np(i)*log(10.))
	    if(count .eq. 3) linefit_sigma(i,count)=sigma_BxStruc(i)/(S2_Bx(i)*log(10.))
	    if(count .eq. 4) linefit_sigma(i,count)=sigma_ByStruc(i)/(S2_By(i)*log(10.))
	    if(count .eq. 5) linefit_sigma(i,count)=sigma_BzStruc(i)/(S2_Bz(i)*log(10.))
	    if(count .eq. 6) linefit_sigma(i,count)=sigma_VxStruc(i)/(S2_Vx(i)*log(10.))
	    if(count .eq. 7) linefit_sigma(i,count)=sigma_VyStruc(i)/(S2_Vy(i)*log(10.))
	    if(count .eq. 8) linefit_sigma(i,count)=sigma_VzStruc(i)/(S2_Vz(i)*log(10.))
            if(count.eq.9)linefit_sigma(i,count)=sigma_BtraceStruc(i)/(Btrace(i)*log(10.))	
            if(count.eq.10)linefit_sigma(i,count)=sigma_VtraceStruc(i)/(Vtrace(i)*log(10.))
	    S(count) = S(count) + (1./(linefit_sigma(i,count)**2))
	    Sx(count) = Sx(count) + taulog(i)/(linefit_sigma(i,count)**2)
	    Sy(count) = Sy(count) + llog(i)/(linefit_sigma(i,count)**2)
	    Sxx(count) = Sxx(count) + (taulog(i)**2)/(linefit_sigma(i,count)**2)
	    Sxy(count) = Sxy(count) + (taulog(i)*llog(i))/(linefit_sigma(i,count)**2)
	  endif
  77	continue
c
	if(isum .lt. 5) then
	  S(count) = bad					
	  Sx(count) = bad
	  Sy(count) = bad
	  Sxx(count) = bad
	  Sxy(count) = bad
	endif
c
c...Compute the fits: y = a + bx
	 if((abs(S(count)-bad) .gt. 1.) .and. (abs(Sx(count)-bad) .gt. 1.) .and.
     .      (abs(Sy(count)-bad) .gt. 1.) .and. (abs(Sxx(count)-bad) .gt. 1.) .and.
     .      (abs(Sxy(count)-bad) .gt. 1.)) then
	 delta_fit(count) = (S(count)*Sxx(count))-(Sx(count)**2)
 	 if(count .eq. 1) then
          Tp_A = ((Sxx(count)*Sy(count)) - (Sx(count)*Sxy(count)))/delta_fit(count)
	    Tp_B = ((S(count)*Sxy(count)) - (Sx(count)*Sy(count)))/delta_fit(count)
          Tp_sigmaA = sqrt(Sxx(count)/delta_fit(count))
	    Tp_sigmaB = sqrt(S(count)/delta_fit(count))
          Tp_sigmaA = (10.0**Tp_A)*log(10.0)*Tp_sigmaA
          Tp_A = 10.**Tp_A
	  endif
	  if(count .eq. 2) then
          Np_A = ((Sxx(count)*Sy(count)) - (Sx(count)*Sxy(count)))/delta_fit(count)
	    Np_B = ((S(count)*Sxy(count)) - (Sx(count)*Sy(count)))/delta_fit(count)
          Np_sigmaA = sqrt(Sxx(count)/delta_fit(count))
	    Np_sigmaB = sqrt(S(count)/delta_fit(count))
          Np_sigmaA = (10.0**Np_A)*log(10.0)*Np_sigmaA
	    Np_A = 10.**Np_A
	  endif
	  if(count .eq. 3) then
          Bx_A = ((Sxx(count)*Sy(count)) - (Sx(count)*Sxy(count)))/delta_fit(count)
	    Bx_B = ((S(count)*Sxy(count)) - (Sx(count)*Sy(count)))/delta_fit(count)
          Bx_sigmaA = sqrt(Sxx(count)/delta_fit(count))
	    Bx_sigmaB = sqrt(S(count)/delta_fit(count))
          Bx_sigmaA = (10.0**Bx_A)*log(10.0)*Bx_sigmaA
	    Bx_A = 10.**Bx_A
	  endif
	  if(count .eq. 4) then
          By_A = ((Sxx(count)*Sy(count)) - (Sx(count)*Sxy(count)))/delta_fit(count)
	    By_B = ((S(count)*Sxy(count)) - (Sx(count)*Sy(count)))/delta_fit(count)
	    By_sigmaA = sqrt(Sxx(count)/delta_fit(count))
	    By_sigmaB = sqrt(S(count)/delta_fit(count))	
 	    By_sigmaA = (10.0**By_A)*log(10.0)*By_sigmaA
          By_A = 10.**By_A
	  endif
	  if(count .eq. 5) then
          Bz_A = ((Sxx(count)*Sy(count)) - (Sx(count)*Sxy(count)))/delta_fit(count)
          Bz_B = ((S(count)*Sxy(count)) - (Sx(count)*Sy(count)))/delta_fit(count)
          Bz_sigmaA = sqrt(Sxx(count)/delta_fit(count)) 
          Bz_sigmaB = sqrt(S(count)/delta_fit(count))
          Bz_sigmaA = (10.0**Bz_A)*log(10.0)*Bz_sigmaA
          Bz_A = 10.**Bz_A
	  endif
	  if(count .eq. 6) then
          Vx_A = ((Sxx(count)*Sy(count)) - (Sx(count)*Sxy(count)))/delta_fit(count)
          Vx_B = ((S(count)*Sxy(count)) - (Sx(count)*Sy(count)))/delta_fit(count)
          Vx_sigmaA = sqrt(Sxx(count)/delta_fit(count))
          Vx_sigmaB = sqrt(S(count)/delta_fit(count))	
          Vx_sigmaA = (10.0**Vx_A)*log(10.0)*Vx_sigmaA
	    Vx_A = 10.**Vx_A
	  endif
	  if(count .eq. 7) then
          Vy_A = ((Sxx(count)*Sy(count)) - (Sx(count)*Sxy(count)))/delta_fit(count)
	    Vy_B = ((S(count)*Sxy(count)) - (Sx(count)*Sy(count)))/delta_fit(count)
          Vy_sigmaA = sqrt(Sxx(count)/delta_fit(count))
	    Vy_sigmaB = sqrt(S(count)/delta_fit(count))
          Vy_sigmaA = (10.0**Vy_A)*log(10.0)*Vy_sigmaA
	    Vy_A = 10.**Vy_A
	  endif
	  if(count .eq. 8) then
          Vz_A = ((Sxx(count)*Sy(count)) - (Sx(count)*Sxy(count)))/delta_fit(count)
	    Vz_B = ((S(count)*Sxy(count)) - (Sx(count)*Sy(count)))/delta_fit(count)
          Vz_sigmaA = sqrt(Sxx(count)/delta_fit(count))
	    Vz_sigmaB = sqrt(S(count)/delta_fit(count))
          Vz_sigmaA = (10.0**Vz_A)*log(10.0)*Vz_sigmaA
	    Vz_A = 10.**Vz_A
	  endif
	  if(count .eq. 9) then
          Btrace_A=((Sxx(count)*Sy(count))-(Sx(count)*Sxy(count)))/delta_fit(count)
	    Btrace_B=((S(count)*Sxy(count))-(Sx(count)*Sy(count)))/delta_fit(count)
	    Btrace_sigmaA = sqrt(Sxx(count)/delta_fit(count))
	    Btrace_sigmaB = sqrt(S(count)/delta_fit(count))	
	    Btrace_sigmaA = (10.0**Btrace_A)*log(10.0)*Btrace_sigmaA
          Btrace_A = 10.**Btrace_A
	 endif
	 if(count .eq. 10) then
          Vtrace_A=((Sxx(count)*Sy(count))-(Sx(count)*Sxy(count)))/delta_fit(count)
	    Vtrace_B=((S(count)*Sxy(count))-(Sx(count)*Sy(count)))/delta_fit(count)
          Vtrace_sigmaA = sqrt(Sxx(count)/delta_fit(count))
	    Vtrace_sigmaB = sqrt(S(count)/delta_fit(count))
          Vtrace_sigmaA = (10.0**Vtrace_A)*log(10.0)*Vtrace_sigmaA
	    Vtrace_A = 10.**Vtrace_A
	 endif
	 else
	 delta_fit(count) = bad
 	 if(count .eq. 1) then
          Tp_A = bad
	    Tp_B = bad
          Tp_sigmaA = bad
	    Tp_sigmaB = bad
	 endif
	 if(count .eq. 2) then
          Np_A = bad
	    Np_B = bad
          Np_sigmaA = bad
	    Np_sigmaB = bad
       endif
	 if(count .eq. 3) then
          Bx_A = bad
	    Bx_B = bad
          Bx_sigmaA = bad
	    Bx_sigmaB = bad
       endif
	 if(count .eq. 4) then
	    By_A = bad
	    By_B = bad
	    By_sigmaA = bad
	    By_sigmaB = bad	
	  endif
	  if(count .eq. 5) then
          Bz_A = bad
	    Bz_B = bad
          Bz_sigmaA = bad
          Bz_sigmaB = bad
	  endif
	  if(count .eq. 6) then
          Vx_A = bad
	    Vx_B = bad
          Vx_sigmaA = bad 
	    Vx_sigmaB = bad	
	  endif
	  if(count .eq. 7) then
          Vy_A = bad
	    Vy_B = bad
          Vy_sigmaA = bad
	    Vy_sigmaB = bad	
	  endif
	  if(count .eq. 8) then
          Vz_A = bad
	    Vz_B = bad
          Vz_sigmaA = bad
	    Vz_sigmaB = bad
	  endif
	  if(count .eq. 9) then
          Btrace_A = bad
	    Btrace_B = bad
	    Btrace_sigmaA = bad
	    Btrace_sigmaB = bad	
	 endif
	 if(count .eq. 10) then
          Vtrace_A = bad
	    Vtrace_B = bad
          Vtrace_sigmaA = bad
	    Vtrace_sigmaB = bad	
	 endif
	endif
c
 1000	continue
	return
	end
	
c... detrend.for ...................................................
c
c  Modified version of below to perform linear fit and detrend of
c    data in Tessein hourly analysis.
c
c  C.W.S. -- 11/2/2007.
c...................................................................
c
c...Linfit.for
c
c Modified to compute linear fits to log/log scaling of the
c   observed dependence of P_perp/P_par upon proton beta.
c
c Uses least-squares routine pulled from Bevington.
c All data equally weighted, although this puts greater emphasis
c on the higher frequencies.
c...................................................................
c
c
c
c
c
c
	subroutine detrend(npts)
c
c.. input data:
	real*8 Bx(60), By(60), Bz(60), Vx(60), Vy(60), Vz(60), Np(60), Tp(60)
	common /SWEPAM/ Bx, By, Bz, Vx, Vy, Vz
c
	integer*4 iyear, icnt_mag(60)
	real*8 time, alpha(60), Vp(60), Vr(60), Vt(60), Vn(60)
	real*8 zlambda(60), delta_dat(60), Br(60), Bt(60), Bn(60), Bmag(60)
	real*8 Xsc(60), Ysc(60), Zsc(60), times
	common /hourly/ time, Np, Tp, alpha, Vp, Vr, Vt, Vn,
     .		Br, Bt, Bn, Bmag, zlambda, delta_dat, Xsc, Ysc, Zsc, 
     .		icnt_mag, times, iyear
c
c.. native variables:
	integer i, npts, idetrend_flag(8)
	real*4 bad
	real*4 delta, a, b, xav, xsqrd, yav, ysqrd, xyav, znpts
	real*4 sigmaa, sigmab, variance
	real*4 ytemp(60), x(60)
	real*4 y(60,8)
c
	bad = -9999.9
c
	do 10 i=1,npts
	  x(i) = float(i)
	  y(i,1) = Bx(i)
	  y(i,2) = By(i)
	  y(i,3) = Bz(i)
	  y(i,4) = Vx(i)
	  y(i,5) = Vy(i)
	  y(i,6) = Vz(i)
	  y(i,7) = Np(i)
	  y(i,8) = Tp(i)
   10	continue
c
	do 500 j=1,8
c
	do 50 i=1,npts
	  ytemp(i) = y(i,j)
   50	continue
c
	isum = 0
	yav = 0.
	xav = 0.
	xyav = 0.
	xsqrd = 0.
	ysqrd = 0.
	do 80 i=1,npts
	  if(abs(ytemp(i)-bad) .gt. 1.) then
		yav = YAV + ytemp(i)
		ysqrd = ysqrd + ytemp(i)**2
		xav = XAV + x(i)
		xsqrd = XSQRD + x(i)**2
		xyav = XYAV + ytemp(i)*x(i)
		isum = isum + 1
	  endif
   80	continue
c
	idetrend_flag(j) = 0
	if(isum .lt. 2) then			! if there is insuficient point
	    idetrend_flag(j) = 1		! to compute a linear fit, then
	    go to 500				! trap out short data sets
	endif
c
c ...and compute statistics of the fit line
c
	znpts = float(isum)
	DELTA = znpts*XSQRD-(XAV)**2
	A = (XSQRD*YAV-XAV*XYAV) / DELTA
	B = (znpts*XYAV-XAV*YAV) / DELTA
	variance = (ysqrd + A*A*znpts + B*B*xsqrd 
     .	   - 2.*(a*yav + B*xyav -A*B*xav)) / (znpts-2.)
	sigmaa = sqrt(variance*xsqrd/delta)
	sigmab = sqrt(variance*znpts/delta)
c
	do 100 i=1,npts
	  if(abs(ytemp(i)-bad) .gt. 1.) then
		yfit = A + B*x(i)
		ytemp(i) = ytemp(i) - yfit
	  endif
	  y(i,j) = ytemp(i)
  100	continue
c
	do 200 i=1,npts
	  Bx(i) = y(i,1) 
	  By(i) = y(i,2)
	  Bz(i) = y(i,3) 
	  Vx(i) = y(i,4) 
	  Vy(i) = y(i,5) 
	  Vz(i) = y(i,6) 
	  Np(i) = y(i,7) 
	  Tp(i) = y(i,8) 
  200	continue
c
  500	continue				! loop for data series
c
	return
	END
