pro aia_mag_merit,fileset,hmax,mag,ds,ff,fov
;+
; Project     : AIA/SDO
;
; Name        : AIA_MAG_LOOP
;
; Category    : SDO/AIA data modeling 
;
; Explanation : Figure of merits
;
; Syntax      : IDL>aia_mag_loop,fileset,hmax,mag,fov,ds
;
; Inputs      : fileset	 = filename 
;		hmax	 = height range of magnetic field model
;
; Outputs:     
;
; History     : 5-May-2011, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;__________________READ LOOPFILE_________________________________
coeff_file=fileset+'_coeff.dat'
readcol,coeff_file,m1,x1,y1,z1
ns      =n_elements(m1)
coeff   =[[m1],[x1],[y1],[z1]]

;__________________READ MAGNETIC FIELD MODEL_____________________
fieldfile =fileset+'*_field_nh.sav'
restore,fieldfile       ;-->field_,curr,dev_,dev_loop,ntwistfit
bnp	=coeff_fit[0]
ntwist	=coeff_fit[1]
lrope	=coeff_fit[2]
polar	=1.

;__________________FIELD-OF-VIEW_________________________________
xc	=(fov[0]+fov[2])/2.
yc	=(fov[1]+fov[3])/2.
nx	=50
ny	=50
nz	=20
x1c	=xc-ds*(nx/2+0.5) 
y1c	=yc-ds*(ny/2+0.5) 
z1c	=1.0
x	=x1c+ds*findgen(nx)
y	=y1c+ds*findgen(ny)
z	=z1c+ds*findgen(nz)
x2c	=x(nx-1)
y2c	=y(ny-1)
z2c	=z(nz-1)
print,'Datacube dimensions = nx,ny,nz,',nx,ny,nz

;___________________MAGNETIC FIELD LINE_________________________
bx_	=fltarr(nx,ny,nz)
by_	=fltarr(nx,ny,nz)
bz_	=fltarr(nx,ny,nz)
b_	=fltarr(nx,ny,nz)
b_tot	=0.
n_tot	=0.
b_min	=0.
b_max	=0.
for iz=0,nz-1 do begin
 for iy=0,ny-1 do begin
  for ix=0,nx-1 do begin
   x1	=x(ix)
   y1	=y(iy)
   z1	=z(iz)
   b_twist=2.*!pi*ntwist/len 
   aia_mag_vector,x1,y1,z1,ds,hmax,coeff,curr,b_twist,polar,x2,y2,z2,b
   bx_(ix,iy,iz)=(x2-x1)/ds
   by_(ix,iy,iz)=(y2-y1)/ds
   bz_(ix,iy,iz)=(z2-z1)/ds
   b_(ix,iy,iz) =b_tot
   b_tot	=b_tot+(abs(b))
   b_max	=b_max > b
   b_min	=b_min < b
  endfor
 endfor
endfor
n_tot	=nx*ny*nz

;_____________________CURRENTS____________________________________ 
jx	=fltarr(nx,ny,nz)
jy	=fltarr(nx,ny,nz)
jz	=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)
   dbx_dz=bx_(ix,iy,iz+1)-bx_(ix,iy,iz-1)
   dby_dx=by_(ix+1,iy,iz)-by_(ix-1,iy,iz)
   dby_dz=by_(ix,iy,iz+1)-by_(ix,iy,iz-1)
   dbz_dx=bz_(ix+1,iy,iz)-bz_(ix-1,iy,iz)
   dbz_dy=bz_(ix,iy+1,iz)-bz_(ix,iy-1,iz)
   jx(ix,iy,iz)=dbz_dy-dby_dz
   jy(ix,iy,iz)=dbx_dz-dbz_dx
   jz(ix,iy,iz)=dby_dx-dbx_dy
  endfor
 endfor
endfor

;_____________________CURL B = j x B_____________________________ 
jxb1=fltarr(nx,ny,nz)
jxb2=fltarr(nx,ny,nz)
jxb3=fltarr(nx,ny,nz)
jxb =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)
   jxb(ix,iy,iz)=jxb1(ix,iy,iz)^2+jxb2(ix,iy,iz)^2+jxb3(ix,iy,iz)^2
  endfor
 endfor
endfor
abs_jxb=sqrt(jxb1(1:nx-2,1:ny-2,1:nz-2)^2+jxb2(1:nx-2,1:ny-2,1:nz-2)^2+jxb3(1:nx-2,1:ny-2,1:nz-2)^2)
abs_j  =sqrt(jx(1:nx-2,1:ny-2,1:nz-2)^2+jy(1:nx-2,1:ny-2,1:nz-2)^2+jz(1:nx-2,1:ny-2,1:nz-2)^2)
abs_b  =b_(1:nx-2,1:ny-2,1:nz-2)
s_j=total(abs_jxb/abs_b)/total(abs_j)

;_____________________FIGURE OF MERIT____________________________
b_avg	=b_tot/float(n_tot)
print,' Average magnetic field = ',b_avg*mag 
print,' Maximum magnetic field = ',b_max*mag 
print,' Force freeness         = ',s_j
ff	=[x1c,y1c,x2c,y2c,s_j]
end