pro nlfff_energy_v3,dir,savefile,test
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : NLFFF_ENERGY 
;
; Category    : Magnetic data modeling
;
; Explanation : Calculating magnetic field vectors in a 3D cube
;
; Syntax      : IDL>nlfff_cube,dir,savefile,test
;
; Inputs      : dir,savefile,test
;
; Outputs     ; updates savefile --> x,y,z,bx,by,bz,b,a
;
; History     : 10-Mar-2012, Version 1 reconstructed by Markus J. Aschwanden
;                2-JUL-2019, new analytical solution NLFFF_VECTOR_V3.PRO 
;		 9-JUL-2019, replace ds-->dpix
;		19-JUL-2019, no isotropic correction 
;
; Contact     : aschwanden@lmsal.com
;-

print,'__________________NLFFF_ENERGY____________________________'
restore,dir+savefile    ;--> PARA
fov0    =input.fov0      ;field-of-view
hmin	=input.hmin
hmax	=input.hmax
filename=strmid(savefile,0,15)
x1_sun  =para.x1_sun
x2_sun  =para.x2_sun
y1_sun  =para.y1_sun
y2_sun  =para.y2_sun
qe_rebin=para.qe_rebin
qe_model=para.qe_model
eps	=1.e-6

;__________________DIMENSIONS_____________________________________
restore,dir+savefile    ;--> para,dateobs,noaa,helpos,instr,wave_
filename=strmid(savefile,0,15)
dim     =size(bzmap)
nx      =dim(1)
ny      =dim(2)
fov     =[x1_sun,y1_sun,x2_sun,y2_sun]
dpix    =float(x2_sun-x1_sun)/float(nx)
x       =x1_sun+dpix*findgen(nx)
y       =y1_sun+dpix*findgen(ny)
nz	=long(hmax/dpix)>1

;___________________COEFF_FILE_____________________________________
coeff0	=coeff
coeff0(*,4)=0.

;___________________MAGNETIC 3D CUBE_______________________________
print,'                                [10^30 erg]'
print,'  iz      Bmax         r        Ep     E_free E_np/E_p'
e_h	=fltarr(nz,3)
b2_p	=0.
b2_f	=0.
nh	=20				;EXTRAPOLATION for ih>20; 2014-MAY-18 
rsun_cm	=6.96e10
unit    =(dpix*rsun_cm/1.e10)^3/(8.*!pi)
for iz=0,nh-1 do begin
 rr 	=1.+iz*dpix
 if (rr ge (1.+hmax)) then goto,endloop
 bmax	=0
 xx	=fltarr(nx*ny)
 yy	=fltarr(nx*ny)
 for iy=0,ny-1 do begin
  i1	=nx*iy
  i2	=nx*(iy+1)-1
  xx[i1:i2]=x
  yy[i1:i2]=y[iy]
 endfor
 zz	=sqrt((rr^2-xx^2-yy^2)>eps)
 x_array=[[xx],[yy],[zz]]
 nlfff_vector_v3,coeff ,x_array,b_	
 nlfff_vector_v3,coeff0,x_array,b0_
 bx_	=b_[*,0]
 by_	=b_[*,1]
 bz_	=b_[*,2]
 bx0_	=b0_[*,0]
 by0_	=b0_[*,1]
 bz0_	=b0_[*,2]
 b      =sqrt(bx_^2+by_^2+bz_^2)
 b0     =sqrt(bx0_^2+by0_^2+bz0_^2)
 vector_product_array,b_,b0_,c,c_norm,angle
 b_phi  =b*abs(sin(angle))
 b2free =b_phi^2			;identical B_z component 
 b2_f	=b2_f+total(b2free)	  
 b2_p   =b2_p+total(b0^2)  
 bmax	=bmax > max(b)
 e_h(iz,0)=rr				;altitude
 e_h(iz,1)=b2_p*unit			;potential energy
 e_h(iz,2)=(b2_p+b2_f)*unit		;nonpotential energy 
 q	=0.
 if (e_h(iz,1) gt 0.) then q=e_h(iz,2)/e_h(iz,1)
 epot	=e_h(iz,1)
 efree	=(e_h(iz,2)-e_h(iz,1))
 ep_est	=0.0
 ef_est	=0.0
 if (iz eq 10) then begin
  h1	=rr-1.
  ep1	=epot
  ef1	=efree
 endif
 if (iz eq 19) then begin
  h2	=rr-1.
  ep2	=epot
  ef2	=efree
  nd	=1000
  pp	=-2.0	
  pf	=-2.0	
  dm_	=0.2*(1.+findgen(nd))/float(nd-1)
  qp_	=(1.-(1.+h2/dm_)^pp)/(1.-(1.+h1/dm_)^pp)
  qf_	=(1.-(1.+h2/dm_)^pf)/(1.-(1.+h1/dm_)^pf)
  dmin	=min(abs(qp_-ep2/ep1),imp)
  dmin	=min(abs(qf_-ef2/ef1),imf)
  dp	=dm_(imp)
  df	=dm_(imf)
  ep0	=ep2/(1.-(1.+h2/dp)^pp)
  ef0	=ef2/(1.-(1.+h2/df)^pf)
 endif
 if (iz ge 20) then begin
  h3	=rr-1.
  ep_est =ep0*(1.-(1.+h3/dp)^pp)
  ef_est =ef0*(1.-(1.+h3/df)^pf)
 endif
;........................................................... 
 if ((iz mod 10) eq 0) or (iz eq nx-1) then $
 print,iz,bmax,rr,epot,efree,q,format='(I4,4I10,f10.3)'
 if (iz eq 0) then begin
  b2_free=fltarr(nx,ny)		;save 2D map of free energy at z=0
  for iy=0,ny-1 do begin
   i1	=nx*iy
   i2	=nx*(iy+1)-1
   b2_free(*,iy)=b2free[i1:i2]
  endfor
 endif
endfor
ENDLOOP:
if (nh eq nz) then begin
 e_p	=e_h(nh-1,1)
 e_np	=e_h(nh-1,2)
endif
if (nh lt nz) then begin
 h_	=dpix*(findgen(nz)+0.5)		;altitude
 ep_h	=ep0*(1.-(1.+h_/dp)^pp)
 ef_h   =ef0*(1.-(1.+h_/df)^pf)
 en_h	=ep_h+ef_h
 e_h(*,0)=1.+h_				;altitude
 e_h(*,1)=ep_h*1.e30			;potential energy
 e_h(*,2)=en_h*1.e30			;nonpotential energy 
 e_free	=ef0*1.e30
 e_p	=ep0*1.e30
 e_np	=e_p+e_free
endif
print,iz,bmax,rr,epot,efree,q,format='(I4,4I10,F10.3)'

;_______________________MAGNETIC ENERGY CORRECTION___________________
q_corr	=1./qe_model                  ;missing magnetic energy correction
e_p     =q_corr*e_p
e_np    =q_corr*e_np
e_free  =(e_np-e_p)

;_______________________SAVE ALPHA PARAMETER_________________________
para.qiso_corr=q_corr
save,filename=dir+savefile,input,para,bzmap,bzfull,bzmodel,coeff,$
     field_loop_det,ns_loop_det,wave_loop_det,$
     field_loop,ns_loop,wave_loop,dev_deg,dev_deg2,$
     field_lines,ns_loop,ns_field,e_h,e_p,e_np,b2_free
print,'parameters saved in file = ',savefile
end