pro magn_misalign,field_lines,nf,nseg,ds,coeff,polar,dev,dev_loop,field_vec,qv
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : MAGN_MISALIGN 
;
; Category    : Magnetic data modeling
;
; Explanation : Subroutine in support of MAGN_FIT.PRO,
;		which calcultes the misalignment angle between an observed
;		coronal loop FIELD_LINES and a theoretical NLFFF model.
;               The fitting constraints can be stereoscopic
;               loop coordinates (qv=0.0), photospheric 3D magnetic
;               field vector data (qv=1.0) or a combination of both.
;
; Syntax      : IDL>magn_misalign,field_lines,nf,nseg,ds,coeff,polar,$
;                                 dev,dev_loop,field_vec,qv
;
; Inputs      : field_lines[NS,3,NF] = 3D coordinates of loop fieldlines
;		nf	= number of field lines
;		nseg	= number of loop segments with misalignment angles
;		ds	= spatial resolution along field line
;		coeff	= magnetic model coefficients [b1,x1,y1,z1,a1]
;		polar	= magnetic polarity (+1,-1) or field direction
;		field_vec=[2,3,NV] = 3D magnetic field vectors at phosotpher
;		qv	= weighting of coronal (qv=0.0) and photospheric 3D
;			  field vectors (qv=1.0)  
;
; Outputs     ; dev	= average misalignment angle of all loops 
;		dev_loop[nf] = average misalignment angles per loop 
;
; History     : 16-May-2011, Version 1 written by Markus J. Aschwanden
;               30-Jan-2012, Version 2 released to SSW
;
; Contact     : aschwanden@lmsal.com
;-

dev_loop=fltarr(nf)           			;misalign angle for each loop
for k=0,nf-1 do begin         			;number of loops
 ind  =where(field_lines[*,2,k] ne 0,ns) 
 xl   =field_lines[0:ns-1,0,k]
 yl   =field_lines[0:ns-1,1,k]
 zl   =field_lines[0:ns-1,2,k]
 dev_  =fltarr(nseg)
 for istep=0,nseg-1 do begin
  im    =long((ns-1)*(1+istep)/float(nseg+1))<(ns-2)
  x1    =xl[im]
  y1    =yl[im]
  z1    =zl[im]
  magn_vector,x1,y1,z1,ds,coeff,polar,x2,y2,z2,bx,by,bz,b,endv
  v1    =[xl[im+1]-xl[im],yl[im+1]-yl[im],zl[im+1]-zl[im]]
  v2    =[x2-x1,y2-y1,z2-z1]
  vector_product,v1,v2,v3,v3_norm,angle_rad
  dev_rad =angle_rad < (!pi-angle_rad)
  dev_(istep)=(180./!pi)*dev_rad      		;misalign angle segm
 endfor
 dev_loop(k)=sqrt(total(dev_^2)/float(nseg))  	;leeast square 
endfor                		      		;for k=0,nf-1
dev1	=sqrt(total(dev_loop^2)/float(nf))      ;unweighted least square
dev2	=0.

;3D vector field fit
if (qv gt 0.) then begin
 dim	=size(field_vec)
 nv	=dim[3]
 dev_vect=fltarr(nv)
 for k=0,nv-1 do begin
  xv	=field_vec[*,0,k]
  yv	=field_vec[*,1,k]
  zv	=field_vec[*,2,k]
  magn_vector,xv[0],yv[0],zv[0],ds,coeff,polar,x2,y2,z2,bx,by,bz,b,endv
  v1    =[xv[1]-xv[0],yv[1]-yv[0],zv[1]-zv[0]]
  v2    =[x2-xv[0],y2-yv[0],z2-zv[0]]
  vector_product,v1,v2,v3,v3_norm,angle_rad
  dev_rad =angle_rad < (!pi-angle_rad)
  dev_vect(k)=(180./!pi)*dev_rad      		;misalign angle vector 
 endfor                		      		;for k=0,nf-1
 dev2	=sqrt(total(dev_vect^2)/float(nv))     	;least square
endif

dev	=(1.-qv)*dev1 + qv*dev2			;relative weighting
end