pro magn_merit,runfile,para
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : MAGN_MERIT 
;
; Category    : Magnetic data modeling
;
; Explanation : Calculates figure of merits of a NLFFF solution:
;		divergence-freeness and force-freeness:	MERIT=[ld,lf,s_j]
;
; Syntax      : IDL>magn_merit,runfile,para
;
; Inputs      : runfile   = input file *_cube.sav, *_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 *_fit.sav
;			    -->bzmap,x,y,field_lines,q_scale,angle,cpu,merit,iter
;
; History     : 28-Nov-2011, Version 1 written by Markus J. Aschwanden
;               30-Jan-2012, Version 2 released to SSW
;
; Contact     : aschwanden@lmsal.com
;-

print,'______________MAGN_MERIT___________________'
dpix	=para[9]

savefile=runfile+'_cube.sav'
restore,savefile        ;-->x,y,z,bx,by,bz,b
nx      =n_elements(x)
ny      =n_elements(y)
nz      =n_elements(z)

;_____________________CURRENTS_____________________ 
jx	=fltarr(nx,ny,nz)
jy	=fltarr(nx,ny,nz)
jz	=fltarr(nx,ny,nz)
div_b	=fltarr(nx,ny,nz)
for iz=1,nz-2 do begin
 for iy=1,ny-2 do begin
  for ix=1,nx-2 do begin
   dbx_dy=(bx(ix,iy+1,iz)-bx(ix,iy-1,iz))/2.
   dbx_dz=(bx(ix,iy,iz+1)-bx(ix,iy,iz-1))/2.
   dby_dx=(by(ix+1,iy,iz)-by(ix-1,iy,iz))/2.
   dby_dz=(by(ix,iy,iz+1)-by(ix,iy,iz-1))/2.
   dbz_dx=(bz(ix+1,iy,iz)-bz(ix-1,iy,iz))/2.
   dbz_dy=(bz(ix,iy+1,iz)-bz(ix,iy-1,iz))/2.
   jx(ix,iy,iz)=(dbz_dy-dby_dz)
   jy(ix,iy,iz)=(dbx_dz-dbz_dx)
   jz(ix,iy,iz)=(dby_dx-dbx_dy)
   dbx_dx=(bx(ix+1,iy,iz)-bx(ix-1,iy,iz))/2.
   dby_dy=(by(ix,iy+1,iz)-by(ix,iy-1,iz))/2.
   dbz_dz=(bz(ix,iy,iz+1)-bz(ix,iy,iz-1))/2.
   div_b(ix,iy,iz)=dbx_dx+dby_dy+dbz_dz
  endfor
 endfor
endfor

;_____________________CURL B = j x B_____________________________ 
jxb1=fltarr(nx,ny,nz)
jxb2=fltarr(nx,ny,nz)
jxb3=fltarr(nx,ny,nz)
for iz=1,nz-2 do begin
 for iy=1,ny-2 do begin
  for ix=1,nx-2 do begin
   jxb1(ix,iy,iz)=jy(ix,iy,iz)*bz(ix,iy,iz)-jz(ix,iy,iz)*by(ix,iy,iz)
   jxb2(ix,iy,iz)=jz(ix,iy,iz)*bx(ix,iy,iz)-jx(ix,iy,iz)*bz(ix,iy,iz)
   jxb3(ix,iy,iz)=jx(ix,iy,iz)*by(ix,iy,iz)-jy(ix,iy,iz)*bx(ix,iy,iz)
  endfor
 endfor
endfor
abs_jxb =sqrt(jxb1^2+jxb2^2+jxb3^2)
abs_j   =sqrt(jx^2+jy^2+jz^2)
s_j	=total(abs_jxb(1:nx-2,1:ny-2,1:nz-2)/b(1:nx-2,1:ny-2,1:nz-2))/$
         total(abs_j(1:nx-2,1:ny-2,1:nz-2))

;____________________Force-freeness______________________________
vol	=float(nx-2)*float(ny-2)*float(nz-2)
lf	=total(abs_jxb(1:nx-2,1:ny-2,1:nz-2)^2/b(1:nx-2,1:ny-2,1:nz-2)^4)/vol

;____________________Divergence-freeness_________________________
ld	=total(div_b(1:nx-2,1:ny-2,1:nz-2)^2/b(1:nx-2,1:ny-2,1:nz-2)^2)/vol

;_____________________FIGURE OF MERIT____________________________
print,' Maximum magnetic field = ',max(abs(b)) 
print,' Divergence freeness    = ',ld
print,' Force freeness         = ',lf
print,' Weighted current       = ',s_j

;____________________SAVE MERIT PARAMETER__________________________
savefile2=runfile+'_fit.sav'
restore,savefile2   ;-->bzmap,x,y,field_lines,q_scale,angle,cpu,merit,iter,dev_loop
merit	=[ld,lf,s_j]
save,filename=savefile2,bzmap,x,y,field_lines,q_scale,angle,cpu,merit,iter,dev_loop
print,'Field lines saved in file = ',savefile2
end