function ring_brightness,r_abs,rho_albedo_phase_index_ring_versus_r_sat,$
		cos_inc_angle,cos_emm_angle, 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
	;NOTE :cos_inc_angle,cos_emm_angle, should be positive

	if_lookup_table=1
	ss_phase_index=3.092

		tau=interpol(rho_albedo_phase_index_ring_versus_r_sat(*,1),$
		rho_albedo_phase_index_ring_versus_r_sat(*,0),r_abs)
		albedo=interpol(rho_albedo_phase_index_ring_versus_r_sat(*,2)$
			,rho_albedo_phase_index_ring_versus_r_sat(*,0),r_abs)

	if (if_lookup_table) then begin
		n_tau=where((tau_lookup-tau)^2 eq min((tau_lookup-tau)^2))

		n_albedo=where((albedo_lookup-albedo)^2 eq min((albedo_lookup-albedo)^2))

		solar_elev=asin(cos_inc_angle)/!pi*180.
		if(if_dark_side) then solar_elev=solar_elev*(-1)
		n_solar_elev=where((solar_elev_lookup-solar_elev)^2 eq min((solar_elev_lookup-solar_elev)^2))

		obs_elev=asin(cos_emm_angle)/!pi*180.
		n_obs_elev=where((obs_elev_lookup-obs_elev)^2 eq min((obs_elev_lookup-obs_elev)^2))

		obs_azim=azimuth_inc_em_ring_plane_deg
		n_obs_azim=where((obs_azim_lookup-obs_azim)^2 eq min((obs_azim_lookup-obs_azim)^2))

		n_ss_phase_index=where((ss_phase_index_lookup-ss_phase_index)^2 eq min((ss_phase_index_lookup-ss_phase_index)^2))

		;stop
		return, ring_phase_lookup_table(n_tau,n_albedo,n_solar_elev,n_obs_elev,n_obs_azim,n_ss_phase_index,1)

	endif else begin
		power_pf_fun=interpol(rho_albedo_phase_index_ring_versus_r_sat(*,3)$
			,rho_albedo_phase_index_ring_versus_r_sat(*,0),r_abs)

			;was before: a formula P(s) = C_n s**power_ph_fun from
			;Luke Dones where s = pi - the phase angle in radians
			;now the one from Dones et. 93
			;the normalization coefficient C_n corresponds to n according to the following table
			; (Dones et al 93)
		c_n_versus_n=[[1.000,1.257,1.860,2.624,3.092,3.301,3.642 ,4.300 ,5.108 ,7.518],$
					  [0.637,0.539,0.375,0.218,0.153,0.130,0.0993,0.0578,0.0289,0.00326]]
		c_n=interpol(c_n_versus_n(*,1),c_n_versus_n(*,0),power_pf_fun)
		;help,c_n,tau
		if (if_dark_side(0) ne 0) then begin
			 mu_0_minus_mu=(cos_emm_angle-cos_inc_angle)
			 if (abs(mu_0_minus_mu) le 1.e-10) then mu_0_minus_mu=1.e-10
			 brightness=albedo/4.*c_n*phase_angle^power_pf_fun*$
			 	(exp(-tau/cos_emm_angle)-exp(-tau/cos_inc_angle))*cos_inc_angle/mu_0_minus_mu
			;	brightness=0.
		endif else begin	;illuminated side
			brightness=albedo/4.*cos_inc_angle/(cos_emm_angle+cos_inc_angle)$
				*(1.-exp(-tau*(1./cos_emm_angle+1./cos_inc_angle)))*c_n*phase_angle^power_pf_fun
		endelse
		;stop
		return, brightness
	endelse
end