pro nlfff_volume,dir,filename,vers,para,run_nr,iev,it,b_cube,e_area
;+
; 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_CUBE_____________________________'
thresh_b=para[30]
hmax	=para[12]
savefile=dir+filename+'_field_'+vers+'.sav'
restore,savefile        ;-->bzmap,field_lines,alpha,q_scale,angle,cpu
dim     =size(bzmap)
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)
nz      =long(hmax/dpix)
z	=1.+dpix*findgen(nz)
print,'Datacube dimensions = nx,ny,nz,',nx,ny,nz

coeff_file=dir+filename+'_coeff_'+vers+'.dat'
readcol,coeff_file,m1,x1,y1,z1,a1,/silent
ncoeff	=n_elements(m1)
a0	=fltarr(ncoeff)
coeff   =[[m1],[x1],[y1],[z1],[a1]]
coeff0   =[[m1],[x1],[y1],[z1],[a0]]

;___________________MAGNETIC 3D CUBE_______________________________
bx      =fltarr(nx,ny,nz)
by      =fltarr(nx,ny,nz)
bz      =fltarr(nx,ny,nz)
b       =fltarr(nx,ny,nz)
for iz=0,nz-1 do begin
 bmax	=0
 for iy=0,ny-1 do begin
  xx    =x
  yy    =x*0.+y[iy]
  zz	=x*0.+z[iz]
  x_array=[[xx],[yy],[zz]]
  nlfff_vector,coeff,x_array,b_
  bx(*,iy,iz)=b_[*,0]
  by(*,iy,iz)=b_[*,1]
  bz(*,iy,iz)=b_[*,2]
  bx_=b_[*,0]
  by_=b_[*,1]
  bz_=b_[*,2]
  b_tot=sqrt(bx_^2+by_^2+bz_^2)
  b(*,iy,iz)=b_tot
  bmax=bmax > max(b_tot)
 endfor
 if (iz mod 10) eq 0 then print,'Magnetic field 3D cube iz/nz=',iz,nz,' bmax=',bmax
 if (bmax lt thresh_b) then goto,end_cube
endfor
END_CUBE: 

;___________________ENERGY DISSIPATION__________________________
if (it eq 0) then b_cube=b
if (it eq 0) then e_area=fltarr(nx,ny)

if (it ge 1) then begin
 thresh_e=thresh_b^2
 b_old  =b_cube
 e_old	=b_old^2
 e_new	=b^2
 de	=(e_new-e_old) > 0.		;only dissipation
 b_cube =b

;3D fractal
 ind	=where(de ge thresh_e,vol) 
 l3	=(nx*ny*nz)^(1./3.)
 d3	=alog10(vol>1)/alog10(l3)

;2D fractal
 area_e	=de[*,*,0]
 for iz=1,nz-1 do area_e=area_e > de[*,*,iz] 
 ind	=where(area_e ge thresh_e,area) 
 l2	=sqrt(nx*ny)
 d2	=alog10(area>1)/alog10(l2)

;2D cumulative area
 e_area =e_area > area_e
 ind	=where(e_area ge thresh_e,area_cum) 
 l2	=sqrt(nx*ny)
 d2cum	=alog10(area_cum>1)/alog10(l2)

 outfile='~/work_global/fractal.txt'
 openw,3,outfile,/append
 printf,3,run_nr,iev,it,l2,l3,d2,d2cum,d3,format='(a,4I6,3f8.3)'
 print   ,run_nr,iev,it,l2,l3,d2,d2cum,d3,format='(a,4I6,4f8.3)'
 close,3
endif
end