pro nlfff_cube_v3,dir,savefile,test,cube
;+
; 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,dir,savefile,test,cube
;
; Inputs      : dir,savefile,test
;               cube=0 (minimal cube with nz=3 layers
;               cube=1 (full cube with nz=zrange/dpix layers
;
; Outputs     ; updates savefile --> x,y,z,bx,by,bz,b,a
;
; History     : 10-Mar-2012, Version 1 reconstructed by Markus J. Aschwanden
;               26-Mar-2015, nz=long(hmax/dpix) < nz_cube
;		 2-Jul-2019, new analytical solution NLFFF_VECTOR_V3.PRO
;
; Contact     : aschwanden@lmsal.com
;-

print,'__________________NLFFF_CUBE_____________________________'
restore,dir+savefile    ;--> PARA
fov0    =input.fov0      ;field-of-view
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
hmax    =0.05        ;maximum altitude (solar radii)
dim     =size(bzmap)
nx      =dim(1)
ny      =dim(2)
z1	=sqrt(1.-x1_sun^2-y1_sun^2)
z2	=sqrt(1.-x1_sun^2-y2_sun^2)
z3	=sqrt(1.-x2_sun^2-y1_sun^2)
z4	=sqrt(1.-x2_sun^2-y2_sun^2)
zmin	=min([z1,z2,z3,z4])
zmax	=max([z1,z2,z3,z4])
print,'zmin,zmax,hmax=',zmin,zmax,hmax
dpix    =float(x2_sun-x1_sun)/float(nx)
x       =x1_sun+dpix*findgen(nx)
y       =y1_sun+dpix*findgen(ny)
zrange	=(zmax-zmin)+hmax
if (cube eq 0) then nz=3
if (cube eq 1) then nz=long(zrange/dpix) 	;full resolution of altitude range
z	=zmin+dpix*findgen(nz)		;z-axis
print,'Datacube dimensions = nx,ny,nz,',nx,ny,nz
print,'z-range=',minmax(z)

;___________________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]
  rr	=sqrt(xx^2+yy^2+zz^2)
  x_array=[[xx],[yy],[zz]]
  nlfff_vector_v3,coeff,x_array,b_
  ind	=where(rr > 1.,nind)
if (nind ge 2) then begin
  bx(ind,iy,iz)=b_[ind,0]
  by(ind,iy,iz)=b_[ind,1]
  bz(ind,iy,iz)=b_[ind,2]
  bx_=b_[ind,0]
  by_=b_[ind,1]
  bz_=b_[ind,2]
  b_tot=sqrt(bx_^2+by_^2+bz_^2)
  b(ind,iy,iz)=b_tot
  bmax=bmax > max(b_tot)
endif
 endfor
 print,'Magnetic field 3D cube iz/nz=',iz,nz,' bmax=',bmax
endfor

;_____________________ALPHA-CUBE___________________________________
a       =fltarr(nx,ny,nz)
for iz=1,nz-2 do begin
 for iy=1,ny-2 do begin
  dbx_dy=(bx(*,iy+1,iz)-bx(*,iy-1,iz))/(2.*dpix)
  dbz_dy=(bz(*,iy+1,iz)-bz(*,iy-1,iz))/(2.*dpix)
  dbx_dz=(bx(*,iy,iz+1)-bx(*,iy,iz-1))/(2.*dpix)
  dby_dz=(by(*,iy,iz+1)-by(*,iy,iz-1))/(2.*dpix)
  dby_dx=(shift(by(*,iy,iz),+1)-shift(by(*,iy,iz),-1))/(2.*dpix)
  dbz_dx=(shift(bz(*,iy,iz),+1)-shift(bz(*,iy,iz),-1))/(2.*dpix)
  jx   =(dbz_dy-dby_dz)
  jy   =(dbx_dz-dbz_dx)
  jz   =(dby_dx-dbx_dy)
  bx_   =bx(*,iy,iz)                                  ;weighted by B^2
  by_   =by(*,iy,iz)
  bz_   =bz(*,iy,iz)
  wx    =bx_^2                                  ;weighted by B^2
  wy    =by_^2
  wz    =bz_^2
  a(*,iy,iz)=total(wx*(jx/bx_)+wy*(jy/by_)+wz*(jz/bz_))/total(wx+wy+wz)
 endfor
endfor

;____________________SAVE DATACUBE____________________________________
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,x,y,z,bx,by,bz,b,a
print,'parameters saved in file = ',savefile
end