pro calculate_ring_brightness_plus_planet,planet,em_angles_rad,field_versus_planet_size,$
				phase_angle_deg,subsolar_azimuth_deg,ring_facing_angle_deg, $
				rho_albedo_phase_index_ring_versus_r_sat, $
				ring_phase_lookup_table,label_lookup,tau_lookup,albedo_lookup,$
				solar_elev_lookup,obs_elev_lookup,obs_azim_lookup,ss_phase_index_lookup,$
				azimuth_inc_em_ring_plane_deg,$
				planet_rings,n_az_ring_integration,n_r_ring_integration,polar_r_versus_equat_r
	;calculates brightness of the rings,
	;puts ring image on top of the crescent image (planet) with all the shadows

	;For the ringshine
	n_az=n_az_ring_integration
	n_r_ring=n_r_ring_integration
	;get the phase angle coefficients for the light reflected off Saturn in the loops below
	@~\lightning\idl\iss\scattered_light\get_barkstorm_ab.pro

	sz=size(planet)
	nx=sz(1)
	ny=sz(2)
	pl_center=[nx/2,ny/2]

;for the shadow on the rings
	@~\lightning\idl\iss\scattered_light\solar_direction.pro
	@~\lightning\idl\iss\scattered_light\norm_ring_plane.pro

;  for ring shadow on the planet
	@~\lightning\idl\iss\scattered_light\get_three_d_planet.pro
	;stop

	nx=sz(1)*field_versus_planet_size
	ny=sz(2)*field_versus_planet_size
	pl_center=[nx/2,ny/2]
	planet_rings=fltarr(nx,ny)
	rings_reflection=planet_rings*0.
	rings_transmission=planet_rings+1.
	planet_large=planet_rings
	cos_em_angles_large=planet_rings
	planet_large(pl_center(0)-sz(1)/2.:pl_center(0)+sz(1)/2.-1 $
		,pl_center(1)-sz(2)/2.:pl_center(1)+sz(2)/2.-1)=planet
	cos_em_angles=cos(em_angles_rad)
	cos_em_angles(where(cos_em_angles eq 1.))=0.
	cos_em_angles_large(pl_center(0)-sz(1)/2.:pl_center(0)+sz(1)/2.-1 $
		,pl_center(1)-sz(2)/2.:pl_center(1)+sz(2)/2.-1)=cos_em_angles
	;tvscl,cos_em_angles_large,2
	;planet_rings=double(planet_rings)

	solar_dot_rings=solar_direction##transpose(norm_ring_plane)
	if_ring_dark_side=(solar_dot_rings(0) le 0.)
	print,'if_ring_dark_side',if_ring_dark_side
	;incidence_rings_deg=acos(solar_dot_rings)/!pi*180
	;print,"incidence angle on rings=",incidence_rings_deg
	cos_inc_angle_rings=abs(solar_dot_rings(0))+1e-7/nx
	abs_sin_ring_facing_angle=abs(sin(ring_facing_angle_deg/180.*!pi)+0.0001/nx)

	r_ring_max=max(rho_albedo_phase_index_ring_versus_r_sat(*,0))
	r_ring_min=min(rho_albedo_phase_index_ring_versus_r_sat(*,0))

	dr=(r_ring_max-r_ring_min)/n_r_ring
	d_lon_rad=360./n_az/180.*2*!pi

;  for if_in_planet_shadow
	pl_radius=1.
	@~\lightning\idl\iss\scattered_light\get_large_to_small_axis_ratio_terminat_plane.pro
	solar_cross_rings=crossp(solar_direction,norm_ring_plane)

	for ix=0,nx-1 do begin
	for iy=0,ny-1 do begin
		;UNITS OF R_Saturn
		r_xy=float([ix+1,iy+1]-pl_center)*2./sz(1)
		r_xy_abs=(r_xy(0)^2+r_xy(1)^2)^.5
		r_xyz=[r_xy(0),r_xy(1),r_xy(1)/(tan(ring_facing_angle_deg/180.*!pi)+0.0000001/nx)]
		;if (norm(r_xyz##transpose(norm_ring_plane)) ge 1.e-5) then print, 'r_xyz##transpose(norm_ring_plane)', r_xyz##transpose(norm_ring_plane)
		r_xyz_abs=(abs(r_xyz##transpose(r_xyz)))^0.5
		;;SHADOWS;;
		if_on_ring=(r_xyz_abs(0) lt r_ring_max) and (r_xyz_abs(0) ge r_ring_min)
		;if (if_on_ring(0)) then print,'r_xyz_abs',r_xyz_abs

		if_on_saturn=(cos_em_angles_large(ix,iy) ne 0.)

		if_behind_planet=(if_on_saturn and (r_xyz(2) lt 0))

		;*******Saturn shadow on rings
		if_saturn_shadow=if_in_planet_shadow(r_xyz,solar_direction,pl_radius, $
			large_to_small_axis_ratio_terminat_plane,solar_cross_rings)

		if (if_on_saturn(0)) then begin

			;******* ring shadow on Saturn

			;r_xyz_planet=[r_xy(0),r_xy(1),(1.-r_xy_abs(0)^2)^0.5]; normalised
			_ix=ix-nx/2.+sz(1)/2.
			_iy=iy-ny/2.+sz(2)/2.
			;print,'_ix=',_ix,'_iy=',_iy,'ix=',ix,'iy=',iy;,'pl_center=',pl_center
			r_xyz_planet=transpose(three_d_planet(_ix,_iy,0:2))
		;	print,'r_xyz_planet',r_xyz_planet
			r_xyz_planet=r_xyz_planet/norm(r_xyz_planet)
			;print,r_xy_abs,'r_xyz_planet=',r_xyz_planet
			r_xyz_plnt_dot_norm_rng=r_xyz_planet##transpose(norm_ring_plane)
			;help,r_xyz_plnt_dot_norm_rng(0)
			proj_ring_plane=r_xyz_planet-solar_direction*r_xyz_plnt_dot_norm_rng(0)/solar_dot_rings(0)
			proj_ring_plane_abs=(proj_ring_plane##transpose(proj_ring_plane))^0.5
			if_upstream=proj_ring_plane##transpose(solar_direction)
			norm_xyz=transpose(three_d_planet(_ix,_iy,3:5))
			ring_shadow=(proj_ring_plane_abs lt r_ring_max) $
				and (proj_ring_plane_abs gt r_ring_min)$
				and (norm_xyz##transpose(solar_direction) ge 0.) $
				and (if_upstream(0) gt 0.)
			;if_on_planet=(r_xy_abs le 1.)
			if (ring_shadow(0)) then begin
				tau=interpol(rho_albedo_phase_index_ring_versus_r_sat(*,1) $
					, rho_albedo_phase_index_ring_versus_r_sat(*,0),proj_ring_plane_abs) $
					/cos_inc_angle_rings
				planet_large(ix,iy)=planet_large(ix,iy)*exp(-tau)
				;print, tau;, proj_ring_plane_abs
			endif

			;******* ring shine on Saturn

			; Each radial interval in rho_ring_versus_r_sat is divided in
			;n_az (see above) azimuthal segments, only segments visible from r_xyz_planet are considered
			;spheric planet; proj_r_planet_ring_plane=r_xyz_planet-r_xyz_plnt_dot_norm_rng(0)*norm_ring_plane
			;r_min_visible=1./norm(proj_r_planet_ring_plane)
			r_xyz_planet=transpose(three_d_planet(_ix,_iy,0:2))/float(sz(1)/2.); units of Saturn R
			r_planet_dot_rings=r_xyz_planet##transpose(norm_ring_plane)
			;print,'r_xyz_planet',r_xyz_planet,'float(sz(0)/2.)',float(sz(1)/2.)
			proj_r_planet_ring_plane=r_xyz_planet-r_planet_dot_rings(0)*norm_ring_plane
			norm_xyz_dot_rings=norm_xyz##transpose(norm_ring_plane)
			norm_xyz_proj_ring_plane=norm_xyz-norm_xyz_dot_rings(0)*norm_ring_plane
			r_min_visible=norm(proj_r_planet_ring_plane)+abs(r_planet_dot_rings(0)*norm_xyz_dot_rings(0)/norm(norm_xyz_proj_ring_plane))
			;print,'r_min_visible',r_min_visible
			if (r_min_visible le r_ring_max $
					;and r_min_visible ge 1. $
					and n_az_ring_integration ne 0. $
					and n_r_ring_integration ne 0.) then begin
				;create a lookup table for az. segment directions
				lon_max=acos(r_min_visible/r_ring_max)
				proj_r_planet_ring_plane_norm=proj_r_planet_ring_plane/norm(proj_r_planet_ring_plane)
				n_lons_total=fix(lon_max/d_lon_rad)+1
				lookup_direction_table=fltarr(2*n_lons_total,3); i_lat,x,y,z of the normalized direction
				proj_perp_r_planet_ring_plane_norm=crossp(proj_r_planet_ring_plane_norm,norm_ring_plane)
				lons=(findgen(2*n_lons_total)-n_lons_total)*d_lon_rad
				for i_lon=0,2*n_lons_total-1 do begin
					lookup_direction_table(i_lon,*)=proj_r_planet_ring_plane_norm*cos(lons(i_lon))+proj_perp_r_planet_ring_plane_norm*sin(lons(i_lon))
				endfor
				;plot,lookup_direction_table(*,0),lookup_direction_table(*,1)
				;wset,1
				;plot,lookup_direction_table(*,0),lookup_direction_table(*,2)
				;wset,0
				;stop
				n_0=fix(n_r_ring*(r_min_visible-r_ring_min)/(r_ring_max-r_ring_min))
				if (r_min_visible le r_ring_min) then n_0=0
				for i_r=n_0, n_r_ring do begin  ;iterate over radial segments
					r_i=r_ring_min+dr*i_r
					lon_max=acos(r_min_visible/r_i)
					n_lons=fix(lon_max/d_lon_rad)+1
					for i_lon=0,2*n_lons-1 do begin  ; iterate over azimutal segments
						r_segment=r_i*transpose(lookup_direction_table(i_lon+(n_lons_total-n_lons),*))
						if (not(if_in_planet_shadow(r_segment,solar_direction,pl_radius, $
									large_to_small_axis_ratio_terminat_plane,solar_cross_rings))) then begin
							r_segment_to_planet=r_segment-r_xyz_planet
							r_segment_to_planet_norm=r_segment_to_planet/norm(r_segment_to_planet)
							cos_emm_angle=r_segment_to_planet_norm##transpose(norm_ring_plane)+1.e-5
							phase_angle=acos(r_segment_to_planet_norm##transpose(solar_direction))
							if_dark_side=(cos_emm_angle(0)*solar_dot_rings(0) ge 0.)
							rings_refl=ring_brightness(r_i(0),rho_albedo_phase_index_ring_versus_r_sat,$
									cos_inc_angle_rings,abs(cos_emm_angle(0)),phase_angle,if_dark_side,$
									ring_phase_lookup_table,label_lookup,tau_lookup,albedo_lookup,$
									solar_elev_lookup,obs_elev_lookup,obs_azim_lookup,ss_phase_index_lookup,$
									azimuth_inc_em_ring_plane_deg)
							d_solid_angle=d_lon_rad*dr/norm(r_segment_to_planet)*abs(cos_emm_angle(0))
							ring_shine=rings_refl*d_solid_angle
							;planet_large(ix,iy)=planet_large(ix,iy)+ring_shine/2/!pi
							incidence_flux=ring_shine
							cos_inc_angle=abs(r_segment_to_planet_norm##transpose(r_xyz_planet))
							cos_emm_angle=cos_em_angles_large(ix,iy)
							phase_angle_rad=acos(-r_segment_to_planet_norm(2));dot-product with [0,0,1] observer direction
							;print,'ring_sine=',ring_shine,'cos_inc_angle',cos_inc_angle,'cos_emm_angle',cos_emm_angle,'phase_angle_rad', phase_angle_rad
							planet_large(ix,iy)=planet_large(ix,iy) $
								+saturn_brightness(incidence_flux,cos_inc_angle(0),cos_emm_angle,phase_angle_rad(0),ab_barkstorm_versus_phase_angle)
						endif
					endfor
				endfor
			endif
		endif
		if ((if_on_ring(0)) and not(if_behind_planet)) then begin
			tau=interpol(rho_albedo_phase_index_ring_versus_r_sat(*,1),rho_albedo_phase_index_ring_versus_r_sat(*,0),r_xyz_abs)$
				/abs_sin_ring_facing_angle
			;print,tau
			;r_xyz=r_xyz/r_xyz_abs(0)
			;print,[ix,iy]
			;print,solar_direction##transpose(r_xyz)
			rings_transmission(ix,iy)=exp(-tau)+1.e-10
			;power_pf_fun=interpol(rho_albedo_phase_index_ring_versus_r_sat(*,3),rho_albedo_phase_index_ring_versus_r_sat(*,0),r_xyz_abs)
			if (not(if_saturn_shadow(0))) then begin
				;if_dark_side=(solar_dot_rings(0) le 0)
				;rings_reflection(ix,iy)=1.
				rings_reflection(ix,iy)=ring_brightness(r_xyz_abs(0),rho_albedo_phase_index_ring_versus_r_sat $
					,abs(cos_inc_angle_rings),abs_sin_ring_facing_angle,(180-phase_angle_deg)/180.*!pi,if_ring_dark_side $
					,ring_phase_lookup_table,label_lookup,tau_lookup,albedo_lookup,$
					solar_elev_lookup,obs_elev_lookup,obs_azim_lookup,ss_phase_index_lookup,$
					azimuth_inc_em_ring_plane_deg)
			endif
		endif
		if (not(if_on_ring(0))) then begin ; to fight a very stange bug - sombrero-like edges
			rings_reflection(ix,iy)=0.
		endif
		if (rings_reflection(ix,iy) le 0) then rings_reflection(ix,iy)=0.
	endfor
	endfor
	;print,'90-em. angle',90-acos(abs_sin_ring_facing_angle)/!pi*180,'90-inc_angle',90-acos(cos_inc_angle_rings)/!pi*180
	tvscl,rings_reflection,1
	;print,'max ,min rings_transmission',max(rings_transmission),min(rings_transmission)
	planet_rings=rings_reflection+planet_large*rings_transmission
	sample=ny/2.;-10
	;plot,acos(cos_em_angles_large(*,sample))/!pi*180,planet_rings(*,sample),ps=1
	plot,planet_rings(*,sample),ps=1
	;print,'max ,min planet_rings',max(planet_rings),min(planet_rings)
	;print,'max ,min rings_reflection',max(rings_reflection),min(rings_reflection)
	;print,'max ,min planet_large',max(planet_large),min(planet_large)

	;help,rings_transmission,planet_rings,rings_reflection,planet_large
end