pro aia_mag_inversion,bz,xp,yp,zp,wp,bm,xm,ym,zm

;calculates model [bm,xm,ym,zm] from observables [bz,xp,yp,zp,wp] 

dm_min	=0.001				;sub-photospheric limit (half MDI pixel)
rho_p	=sqrt(xp^2+yp^2)
alpha	=abs(atan(rho_p/zp))		;approximate value
for iter=0,0 do begin
 if (alpha eq 0) then beta=0.
 if (alpha gt 0) then beta=atan((sqrt(9.+8.*(tan(alpha))^2.)-3.)/(4.*tan(alpha)))
 if (alpha ge !pi/2.) then beta=atan(1./sqrt(2.))
 bm	=bz/(cos(beta)^2*cos(alpha-beta))
 beta2  =acos((cos(beta)^3./2.)^(1./3.))
 dm     =wp/(tan(beta2)*cos(alpha)*(1.-0.1*alpha))
 dm	=dm > dm_min			;NEW Aug 9, 2012 
 rm	=1.-dm
 rho_m	=rho_p-dm*sin(alpha-beta)/cos(beta)
 if (rm ge rho_m) then zm=sqrt(rm^2-rho_m^2)
 if (rm le rho_m) then zm=dm_min	;avoiding singularity APRIL 17,2014
 alpha	=abs(atan(rho_m/zm))		;next iteration
endfor
if (xp eq 0) then gamma=!pi/2.
if (xp ne 0) then gamma=abs(atan(yp/xp))
if (xp eq 0) then xsign=1.
if (yp eq 0) then ysign=1.
if (xp ne 0) then xsign=xp/abs(xp)
if (yp ne 0) then ysign=yp/abs(yp)
xm	=rho_m*cos(gamma)*xsign
ym	=rho_m*sin(gamma)*ysign
end