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 Bvar c Uncertainty in Bvar_sigma c Vvar c Uncertainty in 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