;batch file called from calculate_emission_angles.pro, calculate_incidence_angles.pro

;create a 3d_planet_norm(n_pixels_x(planet),n_pixels_y(planet),6)
;with 3d vectors of half-ellipse and 3d norms at these points in global system

sz=size(planet)
nx=sz(1)
ny=sz(2)

;print,'nx,ny',nx,ny

three_d_planet=fltarr(nx,ny,6)
;help,three_d_planet

;x-east in the imageplane,y-north,z-toward the spacecraft
;@~\lightning\idl\iss\scattered_light\get_x_equat_norm_y_equat_norm.pro
z_equat_norm=norm_ring_plane
;print,'x_equat_norm=',x_equat_norm,'y_equat_norm=',y_equat_norm,'z_equat_norm=',z_equat_norm
;print,'crossp(x_equat_norm,y_equat_norm)=',crossp(x_equat_norm,y_equat_norm)

z_eq_plane_dot_z=z_equat_norm##transpose([0.,0.,1.])
;print,'z_eq_plane_dot_z=',z_eq_plane_dot_z
sin_z_eq_plane_dot_z=-(1.-z_eq_plane_dot_z(0)^2)^0.5

;See the notes on solution for the image-plane-to-ellipsoid projection

epsilon_squared=1./polar_r_versus_equat_r^2

;plot,[0.,nx],[0,ny]-pl_center,ps=1
;plot,[0.,nx],[0,1],ps=1
;print,'pl_center',pl_center

for ix=0,nx-1 do begin
for iy=0,ny-1 do begin
	r_xy=float([ix+1.,iy+1.]-pl_center)/pl_center(0)
	;coeff_A=r_xy(1)*(sin_z_eq_plane_dot_z-z_eq_plane_dot_z(0)^2/sin_z_eq_plane_dot_z)
	;coeff_A=r_xy(1)*(sin_z_eq_plane_dot_z+z_eq_plane_dot_z(0)^2/sin_z_eq_plane_dot_z)
	quadr_eq_a=epsilon_squared*z_eq_plane_dot_z(0)^2 +sin_z_eq_plane_dot_z^2
	quadr_eq_b=-2.*r_xy(1)*z_eq_plane_dot_z(0)*sin_z_eq_plane_dot_z*(epsilon_squared-1)
	quadr_eq_c=r_xy(1)^2*(epsilon_squared*sin_z_eq_plane_dot_z^2+z_eq_plane_dot_z(0)^2) $
		+r_xy(0)^2-1.;pl_center(0)^2
	discriminant=quadr_eq_b(0)^2-4.*quadr_eq_a(0)*quadr_eq_c(0)
	if(discriminant ge 0.) then begin
		;print,',r_xy,r_xy(0)^2,discriminant',r_xy,r_xy(0)^2,discriminant
		;print,'quadr_eq_a',quadr_eq_a,'quadr_eq_b',quadr_eq_b,'quadr_eq_c',quadr_eq_c
		;stop
		z=(quadr_eq_b(0)+discriminant^0.5)/2./quadr_eq_a(0)
		point_on_ellipse=[r_xy(0),r_xy(1),z]
		three_d_planet(ix,iy,0:2)=point_on_ellipse*pl_center(0)
		z_equat_point_on_ellipse=point_on_ellipse##transpose(z_equat_norm)
		norm_planet=point_on_ellipse+ $
		;z_equat_norm*z_equat_point_on_ellipse(0)*(1./polar_r_versus_equat_r-1.)
		z_equat_norm*z_equat_point_on_ellipse(0)*(epsilon_squared-1.)
		norm_planet_abs=norm(norm_planet)
		norm_planet=norm_planet/norm_planet_abs(0)
		three_d_planet(ix,iy,3:5)=norm_planet
		;stop
	endif
endfor
;oplot,findgen(ny),three_d_planet(ix,*,2),color=2*ix
;oplot,findgen(ny),three_d_planet(ix,*,1),color=5*ix
;oplot,findgen(ny),three_d_planet(ix,*,4),color=2*ix
;oplot,findgen(ny),three_d_planet(ix,*,5),color=5*ix
;wait,2
endfor
;stop

;tvscl,three_d_planet(*,*,0);,9
;tvscl,three_d_planet(*,*,1),10
;tvscl,three_d_planet(*,*,2),11
;tvscl,three_d_planet(*,*,3),22
;tvscl,three_d_planet(*,*,4),23
;tvscl,three_d_planet(*,*,5),24
;print,1./polar_r_versus_equat_r
;stop