pro magn_cube,runfile,para
;+
; 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,runfile,para
;
; Inputs      : runfile   = input file *_fit.sav
;               para      = input parameters
;                           [ds,ngrid,nmag,nsm,qloop,nseg,acc,amax,hmax,dpix,
;                           niter,qzone,eps,qv,thresh,q_mag,halt,iopt]
;
; 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,'____________MAGN_CUBE______________________________________'
ds      =para[0]
dpix	=para[9]
hmax	=para[8]

;__________________DIMENSIONS_____________________________________
savefile=runfile+'_fit.sav'
restore,savefile        ;-->bzmap,x,y,field_lines,alpha,q_scale,angle,cpu
nx      =n_elements(x)
ny      =n_elements(y)
nz      =long(hmax/dpix)
z	=1.+dpix*findgen(nz)
print,'Datacube dimensions = nx,ny,nz,',nx,ny,nz
coeff_file=runfile+'_coeff.dat'
readcol,coeff_file,m1,x1,y1,z1,a1
coeff   =[[m1],[x1],[y1],[z1],[a1]]

;___________________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]
   b(*,iy,iz)=b_[*,3]
  bmax=bmax > max(b_[*,3])
 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=savefile,para,fov,x,y,bzmap,field_loops,field_select,$
      angle,dev_loop,iter,cpu,field_lines,z,bx,by,bz,b,a,dev_loop2
print,'parameters saved in file = ',savefile
end