pro nlfff_area,dir,savefile,para,test
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : NLFFF_AREA
;
; Category    : Magnetic data modeling
;
; Explanation : Magnetic data modeling: determines area of free energy
;               and area of energy dissipation (decreases of free energy)
;
; Syntax      : IDL>nlfff_area,...
;
; Inputs      : runfile   = input file *_fit.sav
;               para      = input parameters
;
; Outputs     ; runfile   = output file NLFFF_AREA.TXT 
;
; History     : 2-Aug-2014, Version 1 reconstructed by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;print,'__________________NLFFF_AREA____________________________'
restore,dir+savefile    ;--> para,dateobs,noaa,helpos,instr,wave_
filename=strmid(savefile,0,15)
xc      =para[ 0]               ;FOV in solar radii
yc      =para[ 1]
fovx    =para[ 2]
fovy    =para[ 3]
x1      =xc-fovx/2.
x2      =xc+fovx/2.
y1      =yc-fovy/2.
y2      =yc+fovy/2.
ds      =para[10]     ;spatial resolution along loop
thresh_b=para[16]
thresh_e=thresh_b^2
qiso_corr=(!pi/2)^2     		;isotropic twist axis correction

;print,'__________________PREVIOUS TIMESTEP_____________________'
it	=1
fileprev=filename			;NOT OPERATIONAL
if (filename eq fileprev) then it=0	;first time step
savefile=dir+fileprev+'_field_'+vers+'.sav'
exist	=findfile(savefile,count=count)
if (count eq 0) then begin
 it	=0
 savefile=dir+filename+'_field_'+vers+'.sav'
endif
restore,savefile        
b2_old	=b2_free*qiso_corr
dim	=size(b2_old)
nx_old	=dim[1]
ny_old	=dim[2]
dpix	=(x2-x1)/float(nx_old)		;solar radii/pixel

;print,'__________________PRESENT TIMESTEP______________________'
savefile=dir+filename+'_field_'+vers+'.sav'
restore,savefile        ;parameters saved from NLFFF_ENERGY.PRO
	;-->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
b2_new	=b2_free*qiso_corr
dim	=size(b2_new)
nx_new	=dim[1]
ny_new	=dim[2]
nx	=nx_old < nx_new
ny	=ny_old < ny_new

;_________________________ACCUMULATED FREE ENERGY _________________________
if (it eq 0) then b2_cum=b2_new[0:nx-1,0:ny-1]  
if (it ne 0) then begin
 dim	=size(b2_cum)
 nx_cum	=dim[1]
 ny_cum	=dim[2]
 nx	=nx < nx_cum
 ny	=ny < ny_cum
 b2_cum =b2_new[0:nx-1,0:ny-1] > b2_cum[0:nx-1,0:ny-1]
endif
b2_max	=max(b2_cum)
thresh1	=thresh_e < (b2_max/2.)
ind	=where(b2_cum ge thresh1,area)
l_efree =sqrt(area)*dpix		;length scale in solar radii

;_________________________ACCUMULATED ENERGY DISSIPATION____________________
if (it eq 0) then de_cum =b2_new[0:nx-1,0:ny-1]  
if (it ne 0) then begin
 de	=b2_old(0:nx-1,0:ny-1)-b2_new(0:nx-1,0:ny-1) ;decrease in free energy
 de_cum	=de_cum(0:nx-1,0:ny-1) > (de > 0)	     ;accum energy diss map
endif
de_max	=max(de_cum)
thresh2	=thresh_e < (de_max/2.)
ind    =where(de_cum ge thresh2,area)
l_diss =sqrt(area)*dpix			;length scale in solar radii

;___________________________DISPLAY_______________________________________
if (test eq 3) and (it ne 0) then begin
 clearplot
 loadct,5
 zoom    =4
 window,0,xsize=nx*zoom*2,ysize=ny*zoom
 !p.position=[0,0,0.5,1]
 !x.range=[x1,x2]
 !y.range=[y1,y2]
 !x.style=1
 !y.style=1
 plot,[0,0],[0,0]
 tv,bytscl(rebin(b2_cum,nx*zoom,ny*zoom),min=0,max=b2_max),0,0
 x	=x1+dpix*findgen(nx)
 y	=y1+dpix*findgen(ny)
 level  =thresh1*[0.9,1,1.1]
 contour,b2_cum,x,y,level=level,/overplot
 !noeras=1

 !p.position=[0.5,0,1,1]
 plot,[0,0],[0,0]
 level  =thresh2*[0.9,1,1.1]
 tv,bytscl(rebin(de_cum,nx*zoom,ny*zoom),min=0,max=de_max),nx*zoom,0
 contour,de_cum,x,y,level=level,/overplot
endif

;_______________________SAVE PARAMETERS_________________________
save,filename=savefile,para,z,bzmap,bzmodel,field_loops,ns_loop,ns_field,$
    angle,angle2,dev_deg,dev_deg2,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,wave_loop,b2_cum,de_cum,l_efree,l_diss
rsun=696. ;Mm
print,'parameters saved in file = ',savefile
end