pro nlfff_area2,dir,filename,vers,para,run_nr,iev,it,b_area,l_cum,de_cum,thresh_e
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : NLFFF_CUBE 
;
; Category    : Magnetic data modeling
;
; Explanation : Calculating magnetic field vectors in a 3D cube
;
; Syntax      : IDL>nlfff_cube,imagefile,para
;
; Inputs      : runfile   = input file *_fit.sav
;               para      = input parameters
;                           [ds,ngrid,nmag,nsm,qloop,nseg,acc,amax,hcube,dpix,
;                           niter,qzone,eps,qv,thresh,q_mag,halt
;
; Outputs     ; runfile   = output file *_cube.sav
;		            --> x,y,z,bx,by,bz,b,a
;
; History     : 10-Mar-2012, Version 1 reconstructed by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;print,'__________________NLFFF_AREA____________________________'
thresh_b=para[30]
hmax	=para[12]
savefile=dir+filename+'_field_'+vers+'.sav'
restore,savefile        
	;-->para,x,y,z,bzmap,bzmodel,field_loops,ns_loop,$
   	;-->ns_field,dx_euv,angle,dev_deg,iter,cpu,field_lines,z,bx,by,bz,b,a,$
   	;-->qb_rebin,qb_model,nmag,n_nlfff,ind_nlfff,features
dim	=size(bx)
nx	=dim[1]
ny	=dim[2]
x1      =para[17]
y1      =para[18]
x2      =para[19]
y2      =para[20]
dpix    =float(x2-x1)/float(nx)
x       =x1+dpix*findgen(nx)
y       =y1+dpix*findgen(ny)
coeff_file=dir+filename+'_coeff_'+vers+'.dat'
readcol,coeff_file,m1,x1,y1,z1,a1,/silent
coeff   =[[m1],[x1],[y1],[z1],[a1]]
aia_mag_map,coeff,1.,x,y,bzmodel          ;---> B map

;___________________ENERGY DISSIPATION__________________________
if (it eq 0) then begin
 b_area=fltarr(nx,ny)
 l_area =0.
 l_cum	=0.
 de_cum	=0.
 thresh_e=0.
endif

if (it ge 1) then begin
 b_old  =b_area
 e_old	=b_old^2
 e_new	=bzmodel^2
 de	=(e_new-e_old) > 0.		;only dissipation

;2D fractal
 if (it eq 1) then thresh_e=max(de) 
 if (it gt 1) and (it le 6) then thresh_e=thresh_e < max(de) 
 print,'Threshold = ',sqrt(thresh_e)
 ind	=where(de ge thresh_e,area) 
 l 	=sqrt(nx*ny)
 l_area =0.
 if (area gt 0) then l_area=sqrt(area)

;2D cumulative area
 de_cum =de_cum > de 
 ind	=where(de_cum ge thresh_e,area_cum) 
 l_cum=0.
 if (area_cum gt 0) then l_cum=sqrt(area_cum)
 dim=size(de)
 nx	=dim(1)
 ny	=dim(2)
 zoom	=5
 tv,bytscl(rebin(de,nx*zoom,ny*zoom))
 wait,1

 string0,4,thresh_b,str
 outfile='~/work_global/emag_area_'+str+'.txt'
 openw,3,outfile,/append
 printf,3,run_nr,iev,it,l,l_area,l_cum,l_cum/l,format='(a,5I6,f8.3)'
 print   ,run_nr,iev,it,l,l_area,l_cum,l_cum/l,format='(a,5I6,f8.3)'
 close,3
endif
end