Reflected light phase curves

picaso/atmsetup.py

@@ -111,12 +111,17 @@ def get_profile_3d(self):
  electrons = False
  for g in range(ng):
  for t in range(nt):
- ilat = list(read_3d.coords['lat'].values.astype(np.float32)).index(np.float32(latitude[t]))
- ilon = list(read_3d.coords['lon'].values.astype(np.float32)).index(np.float32(longitude[g]))
+ if 'lat2d' in read_3d.coords: # For phase curve, use lat2d and lon2d
+ ilat = list(read_3d.coords['lat2d'].values.astype(np.float32)).index(np.float32(latitude[t]))
+ ilon = list(read_3d.coords['lon2d'].values.astype(np.float32)).index(np.float32(longitude[g]))
+ else: # For normal 3D calculations, use lat and lon
+ ilat = list(read_3d.coords['lat'].values.astype(np.float32)).index(np.float32(latitude[t]))
+ ilon = list(read_3d.coords['lon'].values.astype(np.float32)).index(np.float32(longitude[g]))
  #read = read_3d[int(latitude[t])][int(longitude[g])].sort_values('pressure').reset_index(drop=True)
  read = read_3d.isel(lon=ilon,lat=ilat).to_pandas().reset_index().drop(['lat','lon'],axis=1).sort_values('pressure')
  if 'phase' in read.keys():
  read=read.drop('phase',axis=1)
+
  #on the first pass look through all the molecules, parse out the electrons and 
  #add warnings for molecules that aren't recognized
  if first:

picaso/disco.py

@@ -34,7 +34,13 @@ def compute_disco(ng, nt, gangle, tangle, phase_angle):
  """
  #this theta is defined from the frame of the downward propagating beam
  cos_theta = cos(phase_angle)
- longitude = arcsin((gangle-(cos_theta-1.0)/(cos_theta+1.0))/(2.0/(cos_theta+1)))
+
+ # This if/else statement allows for full 0-360 phase curve calculations for the reflected case
+ if phase_angle <= pi:
+ longitude = arcsin((gangle-(cos_theta-1.0)/(cos_theta+1.0))/(2.0/(cos_theta+1)))
+ else: 
+ longitude = - arcsin((gangle-(cos_theta-1.0)/(cos_theta+1.0))/(2.0/(cos_theta+1)))
+
  colatitude = arccos(tangle)#colatitude = 90-latitude 
  latitude = pi/2-colatitude
  f = sin(colatitude) #define to eliminate repitition

picaso/fluxes.py

@@ -416,7 +416,7 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
  matrix of cosine of the incident angle from geometric.json
  ubar1 : ndarray of float 
  matrix of cosine of the observer angles
- cos_theta : float 
+ cos_theta: float 
  Cosine of the phase angle of the planet 
  F0PI : array 
  Downward incident solar radiation
@@ -439,7 +439,6 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
  constant_forward : float 
  (Optional), If using the TTHG phase function. Must specify the assymetry of forward scatterer. 
  Remember, the output of A & M code does not separate back and forward scattering.
-
  Returns
  -------
  intensity at the top of the atmosphere for all the different ubar1 and ubar2 
@@ -461,11 +460,13 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
 
  #terms not dependent on incident angle
  sq3 = sqrt(3.)
-
+ 
  #================ START CRAZE LOOP OVER ANGLE #================
+
  for ng in range(numg):
  for nt in range(numt):
-
+ u1 = abs(ubar1[ng,nt]) #absolute value here becuase naegative u values causes overflow errors in x
+ u0 = abs(ubar0[ng,nt]) #absolute value here becuase naegative u values causes overflow errors in x
  #get needed chunk for 3d inputs
  #should consider thinking of a better method for when doing 1d only
  cosb = cosb_3d[:,:,ng,nt]
@@ -490,23 +491,23 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
  g2 = (sq3*w0*0.5)*(1.-ftau_cld*cosb) #table 1
  lamda = sqrt(g1**2 - g2**2) #eqn 21
  gama = (g1-lamda)/g2 #eqn 22
- g3 = 0.5*(1.-sq3*ftau_cld*cosb*ubar0[ng, nt]) #table 1 #ubar is now 100x 10 matrix.. 
+ g3 = 0.5*(1.-sq3*ftau_cld*cosb*u0) #table 1 #ubar is now 100x 10 matrix.. 
 
  # now calculate c_plus and c_minus (equation 23 and 24)
  g4 = 1.0 - g3
- denominator = lamda**2 - 1.0/ubar0[ng, nt]**2.0
+ denominator = lamda**2 - 1.0/u0**2.0
 
  #everything but the exponential 
- a_minus = F0PI*w0* (g4*(g1 + 1.0/ubar0[ng, nt]) +g2*g3 ) / denominator
- a_plus = F0PI*w0*(g3*(g1-1.0/ubar0[ng, nt]) +g2*g4) / denominator
+ a_minus = F0PI*w0* (g4*(g1 + 1.0/u0) +g2*g3 ) / denominator
+ a_plus = F0PI*w0*(g3*(g1-1.0/u0) +g2*g4) / denominator
 
  #add in exponential to get full eqn
  #_up is the terms evaluated at lower optical depths (higher altitudes)
  #_down is terms evaluated at higher optical depths (lower altitudes)
- x = exp(-tau[:-1,:]/ubar0[ng, nt])
+ x = exp(-tau[:-1,:]/u0) 
  c_minus_up = a_minus*x #CMM1
  c_plus_up = a_plus*x #CPM1
- x = exp(-tau[1:,:]/ubar0[ng, nt])
+ x = exp(-tau[1:,:]/u0)
  c_minus_down = a_minus*x #CM
  c_plus_down = a_plus*x #CP
 
@@ -518,10 +519,9 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
  exptrm_positive = exp(exptrm) #EP
  exptrm_minus = 1.0/exptrm_positive#exp(-exptrm) #EM
 
-
  #boundary conditions 
  b_top = 0.0 
- b_surface = 0. + surf_reflect*ubar0[ng, nt]*F0PI*exp(-tau[-1, :]/ubar0[ng, nt])
+ b_surface = 0. + surf_reflect*u0*F0PI*exp(-tau[-1, :]/u0)
 
  #Now we need the terms for the tridiagonal rotated layered method
  #if tridiagonal==0:
@@ -570,13 +570,13 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
  #this is a decent assumption because our second order legendre polynomial 
  #is forced to be equal to the rayleigh phase function
  ubar2 = 0.767 # 
- multi_plus = (1.0+1.5*ftau_cld*cosb*ubar1[ng,nt] #!was 3
- + gcos2*(3.0*ubar2*ubar2*ubar1[ng,nt]*ubar1[ng,nt] - 1.0)/2.0)
- multi_minus = (1.-1.5*ftau_cld*cosb*ubar1[ng,nt] 
- + gcos2*(3.0*ubar2*ubar2*ubar1[ng,nt]*ubar1[ng,nt] - 1.0)/2.0)
+ multi_plus = (1.0+1.5*ftau_cld*cosb*u1 #!was 3
+ + gcos2*(3.0*ubar2*ubar2*u1*u1 - 1.0)/2.0)
+ multi_minus = (1.-1.5*ftau_cld*cosb*u1 
+ + gcos2*(3.0*ubar2*ubar2*u1*u1 - 1.0)/2.0)
  elif multi_phase ==1:#'N=1':
- multi_plus = 1.0+1.5*ftau_cld*cosb*ubar1[ng,nt] 
- multi_minus = 1.-1.5*ftau_cld*cosb*ubar1[ng,nt]
+ multi_plus = 1.0+1.5*ftau_cld*cosb*u1 
+ multi_minus = 1.-1.5*ftau_cld*cosb*u1
  ################################ END OPTIONS FOR MULTIPLE SCATTERING####################
 
  G=w0*positive*(multi_plus+gama*multi_minus)
@@ -591,6 +591,7 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
  #define f (fraction of forward to back scattering), 
  #g_forward (forward asymmetry), g_back (backward asym)
  #needed for everything except the OTHG
+
  if single_phase!=1: 
  g_forward = constant_forward*cosb_og
  g_back = constant_back*cosb_og#-
@@ -614,7 +615,8 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
  #rayleigh phase function
  (gcos2))
  elif single_phase==1:#'OTHG':
- p_single=(1-cosb_og**2)/sqrt((1+cosb_og**2+2*cosb_og*cos_theta)**3) 
+ p_single=(1-cosb_og**2)/sqrt((1+cosb_og**2+2*cosb_og*cos_theta)**3)
+
  elif single_phase==2:#'TTHG':
  #Phase function for single scattering albedo frum Solar beam
  #uses the Two term Henyey-Greenstein function with the additiona rayleigh component 
@@ -624,6 +626,7 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
  #second term of TTHG: backward scattering
  +(1-f)*(1-g_back**2)
  /sqrt((1+g_back**2+2*g_back*cos_theta)**3))
+
  elif single_phase==3:#'TTHG_ray':
  #Phase function for single scattering albedo frum Solar beam
  #uses the Two term Henyey-Greenstein function with the additiona rayleigh component 
@@ -639,22 +642,22 @@ def get_reflected_3d(nlevel, wno,nwno, numg,numt, dtau_3d, tau_3d, w0_3d, cosb_3
  ################################ END OPTIONS FOR DIRECT SCATTERING####################
 
  for i in range(nlayer-1,-1,-1):
- #direct beam
- #note when delta-eddington=off, then tau_single=tau, cosb_single=cosb, w0_single=w0, etc
- xint[i,:] =( xint[i+1,:]*exp(-dtau[i,:]/ubar1[ng,nt])
+ xint[i,:] =( xint[i+1,:]*exp(-dtau[i,:]/u1)
  #single scattering albedo from sun beam (from ubar0 to ubar1)
  +(w0_og[i,:]*F0PI/(4.*pi))*
- (p_single[i,:])*exp(-tau_og[i,:]/ubar0[ng,nt])*
- (1. - exp(-dtau_og[i,:]*(ubar0[ng,nt]+ubar1[ng,nt])/(ubar0[ng,nt]*ubar1[ng,nt])))*
- (ubar0[ng,nt]/(ubar0[ng,nt]+ubar1[ng,nt]))
+ (p_single[i,:])*exp(-tau_og[i,:]/u0)*
+ (1. - exp(-dtau_og[i,:]*(u0+u1)/(u0*u1)))*
+ (u0/(u0+u1))
  #three multiple scattering terms 
- +A[i,:]* (1. - exp(-dtau[i,:] *(ubar0[ng,nt]+1*ubar1[ng,nt])/(ubar0[ng,nt]*ubar1[ng,nt])))*
- (ubar0[ng,nt]/(ubar0[ng,nt]+1*ubar1[ng,nt]))
- +G[i,:]*(exp(exptrm[i,:]*1-dtau[i,:]/ubar1[ng,nt]) - 1.0)/(lamda[i,:]*1*ubar1[ng,nt] - 1.0)
- +H[i,:]*(1. - exp(-exptrm[i,:]*1-dtau[i,:]/ubar1[ng,nt]))/(lamda[i,:]*1*ubar1[ng,nt] + 1.0))
+ +A[i,:]* (1. - exp(-dtau[i,:] *(u0+1*u1)/(u0*u1)))*
+ (u0/(u0+1*u1))
+ +G[i,:]*(exp(exptrm[i,:]*1-dtau[i,:]/u1) - 1.0)/(lamda[i,:]*1*u1 - 1.0)
+ +H[i,:]*(1. - exp(-exptrm[i,:]*1-dtau[i,:]/u1))/(lamda[i,:]*1*u1 + 1.0))
  #thermal
 
- xint_at_top[ng,nt,:] = xint[0,:] 
+ #xint = np.nan_to_num(xint, nan=0)
+ xint_at_top[ng,nt,:] = xint[0,:] 
+
  return xint_at_top
 
 @jit(nopython=True, cache=True)
@@ -783,8 +786,8 @@ def get_reflected_1d_deprecate(nlevel, wno,nwno, numg,numt, dtau, tau, w0, cosb,
  #================ START CRAZE LOOP OVER ANGLE #================
  for ng in range(numg):
  for nt in range(numt):
- u1 = ubar1[ng,nt]
- u0 = ubar0[ng,nt]
+ u1 = abs(ubar1[ng,nt]) #absolute value here becuase naegative u values causes overflow errors in x
+ u0 = abs(ubar0[ng,nt]) #absolute value here becuase naegative u values causes overflow errors in x
  if toon_coefficients == 1 : #eddington
  g3 = (2-3*ftau_cld*cosb*u0)/4#0.5*(1.-sq3*cosb*ubar0[ng, nt]) # #table 1 #ubar has dimensions [gauss angles by tchebyshev angles ]
  elif toon_coefficients == 0 :#quadrature
@@ -802,9 +805,11 @@ def get_reflected_1d_deprecate(nlevel, wno,nwno, numg,numt, dtau, tau, w0, cosb,
  #_up is the terms evaluated at lower optical depths (higher altitudes)
  #_down is terms evaluated at higher optical depths (lower altitudes)
  x = exp(-tau[:-1,:]/u0)
+ #x = exp(-tau[:-1,:]/abs(u0))
  c_minus_up = a_minus*x #CMM1
  c_plus_up = a_plus*x #CPM1
  x = exp(-tau[1:,:]/u0)
+ #x = exp(-tau[1:,:]/abs(u0))
  c_minus_down = a_minus*x #CM
  c_plus_down = a_plus*x #CP
 
@@ -816,7 +821,6 @@ def get_reflected_1d_deprecate(nlevel, wno,nwno, numg,numt, dtau, tau, w0, cosb,
  exptrm_positive = exp(exptrm) #EP
  exptrm_minus = 1.0/exptrm_positive#EM
 
-
  #boundary conditions 
  #b_top = 0.0 
 
@@ -1348,6 +1352,9 @@ def get_reflected_1d(nlevel, wno,nwno, numg,numt, dtau, tau, w0, cosb,gcos2, fta
  0.75*(1+cos_theta**2.0) #rayleigh phase function
  )
  )
+ elif single_phase==4:#'LAB':
+ #Phase function as measured by ExCESS, extrapolated to full viewing angles
+ p_single=(ftau_cld*cos_theta)
  #exploring.... 
  #elif single_phase==4:#'P(HG) exact w/ approx costheta'
  # deltaphi=0
@@ -1366,7 +1373,6 @@ def get_reflected_1d(nlevel, wno,nwno, numg,numt, dtau, tau, w0, cosb,gcos2, fta
  # )
  # )
 
-
  ################################ end options for direct scattering ####################
 
  #in codes like DISORT the single and multiple scattering beams are reported separately