pro nlfff_evol,runs_,fov0_,lon_max,tmargin,dt,dir,efile,catfile,evolfile,iev1,iev2,instr ;+ ; Project : SOHO/MDI, SDO/HMI, STEREO ; ; Name : NLFFF_EVOL ; ; Category : Magnetic data modeling ; ; Explanation : Synthesizing and extracting parameters ; of magnetic field evolution from all runs ; ; Syntax : IDL>nlfff_evol, ; ; Inputs : runs =names of runs to average, e.g., runs=['1','2','3'] ; fov0_ =FOV of runs [solar radii], e.g., [0.35,0.35,0.35] ; lon_max =maximum longitude [deg], e.g., lon_max=45 ; tmargin =time margin of flare time interval, tmargin=0.5 hrs ; dt =time step of time series, dt=0.1 hrs ; efile =GOES event catalog: 'goes_noaa_M1.txt' ; ; Outputs ; catfile =output catalog file: instr+'_evol.dat' ; ; evolfile=output save file: instr+'_evol.sav' ; ; History : 4-AUG-2014, Version 1 reconstructed by Markus J. Aschwanden ; ; Contact : aschwanden@lmsal.com ;- ;________________________NUMBER OF EVENTS________________________ readcol,efile,nr_,tstart_,tpeak_,tend_,goes_,noaa_,pos_,$ format='(I,A,A,A,A,A,A)' nev =iev2-iev1+1 print,'Number of events NEV = ',nev ;________________________OUTPUT ARRAYS___________________________ nrun =n_elements(runs_) nt0 =100 npara =31 nt_ =lonarr(nev) iev_ =lonarr(nev) time_ =strarr(nt0,nev) hours_ =fltarr(nt0,nev) nstruc_ =fltarr(nt0,nev,nrun,8) ;features=fltarr[8] nloop_ =fltarr(nt0,nev,nrun) nlfff_ =fltarr(nt0,nev,nrun) amis_ =fltarr(nt0,nev,nrun) ep_ =fltarr(nt0,nev,nrun) efree_ =fltarr(nt0,nev,nrun) ediss_ =fltarr(nt0,nev,nrun) fdiss_ =fltarr(nt0,nev,nrun) len_ =fltarr(nt0,nev,nrun) magn_ =fltarr(nt0,nev,nrun) para_ =fltarr(npara,nev,nrun) valid =lonarr(nev) ;_________________________READ GOES CATALOG__________________ lon =long(strmid(pos_,1,2)) lat =long(strmid(pos_,4,2)) angle_deg=sqrt(lon^2+lat^2) angle_rad=angle_deg*(!pi/180.) a_proj =cos(angle_rad) ;_________________READ AND EXTRACTS FILES__________________________________ for irun=0,nrun-1 do begin run_nr =runs_[irun] fov0 =fov0_[irun] for iev=iev1,iev2 do begin ii =iev-iev1 nlfff_event,efile,iev,dt,tmargin,nt,dateobs_,fov0,fov_,goes,noaa,pos,ang if (long(strmid(pos,4,2)) gt lon_max) then goto,skip_event valid[ii]=1 nt_(ii)=nt iev_(ii)=iev for it=0,nt-1 do begin filename=dateobs_(it) search=dir+dateobs_(it)+'_field_'+instr+'_'+runs_[irun]+'.sav' datafile=findfile(search,count=count) if (count ge 1) then begin print,'Datafile not found: search=',search features=lonarr(8) restore,datafile(0) print,'READ: ',datafile(0),irun,iev,it ;............................................................................ ; save,filename=savefile,para,z,bzmap,bzmodel,field_loops,ns_loop,ns_field,$ ; angle,dev_deg,cpu,field_lines,z,bx,by,bz,b,a,merit,dx_euv,nmag,$ ; qb_rebin,qb_model,qiso_corr,e_h,e_p,e_np,iter,b2_free,n_nlfff,$ ; ind_nlfff,features,q_thresh,b2_cum,de_cum,l_efree,l_diss ;..............................;.............................................. PARA: npa =n_elements(para) para_(0:npa-1,ii,irun)=para(0:npa-1) for j=0,7 do nstruc_(it,ii,irun,j)=features[j] nloop_(it,ii,irun)=n_elements(ns_loop) nlfff_(it,ii,irun)=n_nlfff amis_(it,ii,irun)=median(dev_deg) ENERGIES: ep_(it,ii,irun)=e_p/1.e30 efree_(it,ii,irun)=(e_np-e_p)/1.e30 TIME: year =strmid(dateobs_(it), 0,4) month=strmid(dateobs_(it), 4,2) day =strmid(dateobs_(it), 6,2) hour =strmid(dateobs_(it), 9,2) min =strmid(dateobs_(it),11,2) sec =strmid(dateobs_(it),13,2) time_(it,ii)=year+'-'+month+'-'+day+'T'+hour+':'+min+':'+sec hours_(it,ii)=float(hour)+float(min)/60.+float(sec)/3600. ;hrs of day LENGTH: rsun=696. ;Mm/solar radius len_(it,ii,irun)=rsun*l_diss ;length scale [Mm] if (it eq 1) then len_(0,ii,irun)=len_(1,ii,irun) ;correct it=0 MAGNETIC_FIELD: magn_(it,ii,irun)=sqrt(max(b2_free)) endif ;volume length scale endfor ;for it=0,nt-1 ENERGY_DISSIPATION: efree =efree_(0:nt-1,ii,irun) for it=1,nt-1 do if (efree(it) eq 0) then efree(it)=efree(it-1) ;correct zero gaps ediss =fltarr(nt) fdiss =fltarr(nt) for it=1,nt-1 do ediss(it)=ediss(it-1)+((efree(it-1)-efree(it))>0.) for it=1,nt-2 do fdiss(it)=ediss(it+1)-ediss(it-1) ediss_(0:nt-1,ii,irun)=ediss fdiss_(0:nt-1,ii,irun)=fdiss SKIP_EVENT: endfor ;for ii=0,nev-1 endfor ;for irun=0,nrun-1 ;_____________________MIDNIGHT CORRECTION_________________________________ for iev=iev1,iev2 do begin ii =iev-iev1 nt =nt_(ii) if (valid[ii] eq 1) then begin hours =hours_[0:nt-1,ii] t0_hr =(anytim(tpeak_[ii] ,/sec) mod 86400)/3600. t1_hr =(anytim(tstart_[ii],/sec) mod 86400)/3600. t2_hr =(anytim(tend_[ii] ,/sec) mod 86400)/3600. if (t2_hr lt t1_hr) then t2_hr=t2_hr+24.0 if (t0_hr lt t1_hr) then t0_hr=t0_hr+24.0 ind =where(hours ne 0) ind_step=where((hours(ind) lt hours(ind(0))),nind_step) if (nind_step ge 1) then begin print,'Before midnight correction : ',hours hours(ind(ind_step))=hours(ind(ind_step))+24.0 print,'After midnight correction : ',hours hours_[0:nt-1,ii]=hours ;updating hours endif endif endfor ;_____________________AVERAGING OVER ALL RUNS_____________________________ ep_avg =fltarr(nt0,nev,2) efree_avg=fltarr(nt0,nev,2) ediss_avg=fltarr(nt0,nev,2) fdiss_avg=fltarr(nt0,nev,2) len_avg =fltarr(nt0,nev,2) magn_avg =fltarr(nt0,nev,2) ene =fltarr(nev,2) peak =fltarr(nev,2) dur =fltarr(nev,2) len =fltarr(nev,2) magn =fltarr(nev,2) for iev=iev1,iev2 do begin ii =iev-iev1 nt =nt_(ii) for it=0,nt-1 do begin statistic,ep_(it,ii,*),e_avg,e_sig ep_avg(it,ii,0)=e_avg ep_avg(it,ii,1)=e_sig/sqrt(nrun-1) ;error of mean value statistic,efree_(it,ii,*),e_avg,e_sig efree_avg(it,ii,0)=e_avg efree_avg(it,ii,1)=e_sig/sqrt(nrun-1) ;error of mean value statistic,ediss_(it,ii,*),e_avg,e_sig ediss_avg(it,ii,0)=e_avg ediss_avg(it,ii,1)=e_sig/sqrt(nrun-1) ;error of mean value statistic,fdiss_(it,ii,*),e_avg,e_sig fdiss_avg(it,ii,0)=e_avg fdiss_avg(it,ii,1)=e_sig/sqrt(nrun-1) ;error of mean value statistic,len_(it,ii,*),l_avg,l_sig len_avg(it,ii,0)=l_avg len_avg(it,ii,1)=l_sig/sqrt(nrun-1) ;error of mean value statistic,magn_(it,ii,*),m_avg,m_sig magn_avg(it,ii,0)=m_avg magn_avg(it,ii,1)=m_sig/sqrt(nrun-1) ;error of mean value endfor ;D,P,E,L values if (nt ge 1) then begin ediss =ediss_avg(0:nt-1,ii,0) sdiss =ediss_avg(0:nt-1,ii,1) frate =fdiss_avg(0:nt-1,ii,0) srate =fdiss_avg(0:nt-1,ii,1) length =len_avg(0:nt-1,ii,0) slength =len_avg(0:nt-1,ii,1) bmagn =magn_avg(0:nt-1,ii,0) nran =10 ene_ran =fltarr(nran) peak_ran =fltarr(nran) len_ran =fltarr(nran) for iran=0,nran-1 do begin seed =iran*100 rho =randomn(seed,nran) ene_ran(iran) =max(ediss+sdiss*rho) peak_ran(iran)=max(frate+srate*rho) len_ran(iran) =max(length+slength*rho) frate_ =frate+srate*rho endfor statistic, ene_ran, ene_avg, ene_sig statistic,peak_ran,peak_avg,peak_sig statistic, len_ran, l_avg, l_sig ene(ii,*) =[ ene_avg, ene_sig] peak(ii,*)=[peak_avg,peak_sig] len(ii,*) =[ l_avg, l_sig] t1_hr =(anytim(tstart_[ii],/sec) mod 86400)/3600. t2_hr =(anytim(tend_[ii] ,/sec) mod 86400)/3600. dh =t2_hr-t1_hr dur(ii,0)=(t2_hr-t1_hr) dur(ii,1)=dt statistic,bmagn,m_avg,m_sig magn(ii,0)=m_avg magn(ii,1)=m_sig endif endfor ;____________________CATALOG_____________________________________________ openw,3,catfile printf,3,' NR EP +/- EFREE +/- DUR +/- ENE +/- PEAK +/- LEN +/- B +/-' print ,' NR EP +/- EFREE +/- DUR +/- ENE +/- PEAK +/- LEN +/- B +/-' for nr_a=iev1,iev2 do begin ii =nr_a-iev1 nt =nt_(ii) ep_a =0. sp_a =0. efree_a=0. sfree_a=0. d =0. e =0. p =0. l =0. m =0. dd =0. ee =0. pp =0. ll =0. mm =0. if (nt ge 1) then begin statistic,ep_avg(0:nt-1,ii,0),ep_a,sp_a statistic,efree_avg(0:nt-1,ii,0),efree_a,sfree_a d =dur(ii,0) e =ene(ii,0) p =peak(ii,0) l =len(ii,0) m =magn(ii,0) dd =dur(ii,1) ee =ene(ii,1) pp =peak(ii,1) ll =len(ii,1) mm =magn(ii,1) endif printf,3,nr_a,ep_a,sp_a,efree_a,sfree_a,d,dd,e,ee,p,pp,l,ll,m,mm,$ format='(I4,4F7.1,2F7.2,8F7.1)' print ,nr_a,ep_a,sp_a,efree_a,sfree_a,d,dd,e,ee,p,pp,l,ll,m,mm,$ format='(I4,4F7.1,2F7.2,8F7.1)' endfor close,3 print,'Catalog file written = ',catfile ;____________________SAVEFILE____________________________________________ save,filename=evolfile,nrun,nev,nt0,npara,para_,time_,nt_,hours_,nstruc_,$ nloop_,nlfff_,amis_,ep_avg,efree_avg,ediss_avg,fdiss_avg,len_avg,magn_avg,$ valid,nr_,tstart_,tpeak_,tend_,goes_,noaa_,pos_,dur,peak,ene,iev_ print,'savefile written = ',evolfile end