diff --git a/docs/notebooks/9g_ReflectedPhaseCurve.ipynb b/docs/notebooks/9g_ReflectedPhaseCurve.ipynb new file mode 100644 index 00000000..53379434 --- /dev/null +++ b/docs/notebooks/9g_ReflectedPhaseCurve.ipynb @@ -0,0 +1,534 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# **Reflected Light Phase Curves**\n", + "\n", + "**Citation: Hamill et al. (2024) : Reflected Light Phase Curves with PICASO: A Kepler-7b Case Study**\n", + "\n", + "Computing a reflected light phase curve is similar to computing 3D reflected light spectra, but now we will compute reflected spectra for an entire grid of phase angles.\n", + "\n", + "From the previous tutorials you should understand:\n", + "\n", + "1. How to convert GCM input to PICASO's required Xarray\n", + "2. How to post-process output to append to GCM ouput\n", + "3. How to run a thermal phase curve and analyze the output\n", + "\n", + "In this tutorial, you will learn:\n", + "\n", + "1. How to compute a reflected light phase curve\n", + "2. How to analyze the output" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "\n", + "from astropy import units as u\n", + "\n", + "import xarray as xr\n", + "\n", + "from picaso import justdoit as jdi\n", + "from picaso import justplotit as jpi\n", + "jpi.output_notebook()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setting Up Reflected Light Phase Curve\n", + "\n", + "Reflected light calculations consider only scattering from the dayside hemisphere. The dayside hemisphere will continually come in or out of view depending on the phase angle.\n", + "\n", + "Let's load in opacity for the visible to near-infrared and the 3D hot Jupiter atmosperic model:\n", + "\n", + "We are resampling our opacities here to make phase curve calculations faster! Use resampling with caution in your own phase curve calculations!!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "opacity = jdi.opannection(wave_range=[0.50,0.90], resample=100) # this is the wavelength band of Kepler" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gcm_out = jdi.HJ_pt_3d(as_xarray=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Add chemistry. \n", + "\n", + "We can use user-input chemistry for this example.\n", + "\n", + "To use post-processed chemical abundances, use the same methods in [Modeling a Phase Curve pt 2 (Robbins-Blanch et al. 2022) : Run Thermal Phase Curve w/ Post-Processed Chemical Abundances.](https://natashabatalha.github.io/picaso/notebooks/9f_PhaseCurves-wChemEq.html)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create coords\n", + "lon = gcm_out.coords.get('lon').values\n", + "lat = gcm_out.coords.get('lat').values\n", + "pres = gcm_out.coords.get('pressure').values\n", + "\n", + "fake_chem_H2O = np.random.rand(len(lon), len(lat),len(pres))*0.1+0.1 # create fake data\n", + "fake_chem_H2 = 1-fake_chem_H2O # create data\n", + "\n", + "# put data into a dataset\n", + "ds_chem = jdi.xr.Dataset(\n", + " data_vars=dict(\n", + " H2O=([\"lon\", \"lat\",\"pressure\"], fake_chem_H2O,{'units': 'v/v'}),\n", + " H2=([\"lon\", \"lat\",\"pressure\"], fake_chem_H2,{'units': 'v/v'}),\n", + "\n", + " ),\n", + " coords=dict(\n", + " lon=([\"lon\"], lon,{'units': 'degrees'}), #required\n", + " lat=([\"lat\"], lat,{'units': 'degrees'}), #required\n", + " pressure=([\"pressure\"], pres,{'units': 'bar'})#required*\n", + " ),\n", + " attrs=dict(description=\"coords with vectors\"),\n", + ")\n", + "all_gcm = gcm_out.update(ds_chem)\n", + "\n", + "all_gcm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "all_gcm['temperature'].isel(pressure=34).plot(x='lon',y='lat')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Set up phase curve grid \n", + "\n", + "Define the number of phases, and the minimum and maximum phase.\n", + "\n", + "For reflected phase curves, the minimum spatial resolution accepted by PICASO is 6 x 6 (num_gangle=6, num_tangle=6)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "case_3d = jdi.inputs()\n", + "\n", + "n_phases = 8\n", + "min_phase = 0\n", + "max_phase = 2*np.pi\n", + "phase_grid = np.linspace(min_phase,max_phase,n_phases)#between 0 - 2pi\n", + "#send params to phase angle routine\n", + "#case_3d.phase_angle(phase_grid=phase_grid,\n", + "# num_gangle=6, num_tangle=6,calculation='reflected')\n", + "case_3d.phase_curve_geometry('reflected', phase_grid=phase_grid, num_gangle=6, num_tangle=6)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A few notes: \n", + "\n", + "- Phase = 0 is the secondary eclipse (maximum reflection)\n", + "- Phase = 2*pi (~6.28) is the primary eclipse (no reflection).\n", + "\n", + "Things to check: \n", + "- Ensure that (0 degrees longitude, 0 degrees latitude) represents your substellar point on the planet. If it does not, you will need to change your reference point in your GCM.\n", + "- For the reflected case, zero_point must be set to 'secondary eclipse'.\n", + "- For a tidally locked planet, set shift to 0 at all phase angles." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "case_3d.atmosphere_4d(all_gcm, shift = np.zeros(n_phases), zero_point='secondary_eclipse',\n", + " plot=True,verbose=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### `Atmosphere_4d`\n", + "\n", + "Atmosphere_4d function works by adding 'phase' as another coordinate in the GCM. 'lat2d' and 'lon2d' are also created and store the specific lats/lons for which the dayside hemisphere is visible to the observer at each individual phase." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "case_3d.inputs['atmosphere']['profile']" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compute the reflected light phase curve.\n", + "\n", + "Note: Planetary radius, orbital distance, and stellar radius is needed to compute reflected phase curves in units of planet/star flux ratio (Fp/Fs)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "case_3d.gravity(radius=1,radius_unit=jdi.u.Unit('R_jup'),\n", + " mass=1, mass_unit=jdi.u.Unit('M_jup')) #any astropy units available\n", + "case_3d.star(opacity,5000,0,4.0, radius=1, radius_unit=jdi.u.Unit('R_sun'), semi_major=0.06, semi_major_unit=u.Unit('au'))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Phase curve computation is below. Computation time will vary widely depending on number of phases and spatial resolution used. We recommend a high performance computer for any phase curves greater than 10x10 spatial resolution. Please see the Appendix of Hamill et al. (2024) for example computation times." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "case_3d.inputs['atmosphere']['profile'].isel(phase=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "allout = case_3d.phase_curve(opacity, n_cpu = 3,#jdi.cpu_count(),\n", + " full_output=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot the phase curve!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "to_plot = 'fpfs_reflected' \n", + "#collapse = [0.40 ]#micron, could select one or more specific wavelengths instead of the mean \n", + "collapse='np.mean'\n", + "\n", + "phases, all_curves, all_ws, fig=jpi.phase_curve(allout, to_plot,\n", + " collapse=collapse, R=10,reorder_output=True)\n", + "\n", + "jpi.show(fig)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Phase-resolved spectra used to create the phase curve.\n", + "\n", + "\n", + "Since chemistry is randomized for this tutorial and there are no clouds, the spectra is mostly symmetrical across 0 degrees phase angle." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#same old same old\n", + "wno =[];fpfs=[];legend=[]\n", + "for iphase in allout.keys():\n", + " w,f = jdi.mean_regrid(allout[iphase]['wavenumber'],\n", + " allout[iphase]['fpfs_reflected'], R=100)\n", + " wno+=[w]\n", + " fpfs+=[f*1e6]\n", + " legend +=[str(int(iphase*180/np.pi))]\n", + "jpi.show(jpi.spectrum(wno, fpfs, plot_width=500,legend=legend,\n", + " palette=jpi.pals.viridis(n_phases)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Reflected phase curve with `Virga` clouds \n", + "\n", + "In reflected light brightness can be heavily influenced by cloud layers. Let's add a user-input cloud layer and compare ouputs**\n", + "\n", + "To use Virga clouds, please use the same steps found in [Modeling a 3D spectrum (Adams et al. 2022) : Run Cloudy 3D spectra](https://natashabatalha.github.io/picaso/notebooks/9c_PostProcess3Dinput-Clouds.html#Post-Process-Clouds:-virga-cloud-model) or [Modeling a Phase Curve pt 2 : Run Thermal Phase Curves w/ Post-Processed Virga Models](https://natashabatalha.github.io/picaso/notebooks/9f_PhaseCurves-wChemEq.html#Run-Thermal-Phase-Curve-w/-Post-Processed-virga-Models)\n", + "\n", + "If you use Virga clouds for reflected phase curves, we highly advise against regridding beforehand. This means your Virga clouds may take a while to compute, but it prevents issues with the re-mapping used in `clouds_4d`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create coords\n", + "lon = gcm_out.coords.get('lon').values\n", + "lat = gcm_out.coords.get('lat').values\n", + "pres = gcm_out.coords.get('pressure').values[:-1]\n", + "wno_grid = np.linspace(33333,1e4,10) #cloud properties are defined on a wavenumber grid\n", + "\n", + "#create box-band cloud model\n", + "fake_opd = np.zeros((len(lon), len(lat),len(pres), len(wno_grid))) # create fake data\n", + "where_lat = np.where(((lat>-50) & (lat<50)))#creating a grey cloud band\n", + "where_lon = np.where(((lon>-90) & (lon<0)))#creating a grey cloud band\n", + "where_pres = np.where(((pres<0.01) & (pres>.00001)))#creating a grey cloud band, 10 mbar to 0.01 mbar\n", + "for il in where_lat[0]:\n", + " for ilon in where_lon[0]:\n", + " for ip in where_pres[0]:\n", + " fake_opd[ilon,il,ip,:]=10 #optical depth of 10 (>>1)\n", + "\n", + "#make up asymmetry and single scattering properties\n", + "fake_asymmetry_g0 = 0.5 + np.zeros((len(lon), len(lat),len(pres), len(wno_grid)))\n", + "fake_ssa_w0 = 0.9 + np.zeros((len(lon), len(lat),len(pres), len(wno_grid)))\n", + "\n", + "# put data into a dataset\n", + "ds_cld= xr.Dataset(\n", + " data_vars=dict(\n", + " opd=([\"lon\", \"lat\",\"pressure\",\"wno\"], fake_opd,{'units': 'depth per layer'}),\n", + " g0=([\"lon\", \"lat\",\"pressure\",\"wno\"], fake_asymmetry_g0,{'units': 'none'}),\n", + " w0=([\"lon\", \"lat\",\"pressure\",\"wno\"], fake_ssa_w0,{'units': 'none'}),\n", + " ),\n", + " coords=dict(\n", + " lon=([\"lon\"], lon,{'units': 'degrees'}),#required\n", + " lat=([\"lat\"], lat,{'units': 'degrees'}),#required\n", + " pressure=([\"pressure\"], pres,{'units': 'bar'}),#required\n", + " wno=([\"wno\"], wno_grid,{'units': 'cm^(-1)'})#required for clouds\n", + " ),\n", + " attrs=dict(description=\"coords with vectors\"),\n", + ")\n", + "\n", + "ds_cld" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sanity Check \n", + "\n", + "Plot of our optically thick cloud layer that lies on the western dayside. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds_cld['opd'].isel(pressure=where_pres[0][0],wno=0).plot(x='lon',y='lat')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Run clouds_4d \n", + "\n", + "Set up phase curve calculation same as before. The only difference is we are now running clouds_4d as well." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "case_3d_clouds = jdi.inputs()\n", + "\n", + "n_phases = 8\n", + "min_phase = 0\n", + "max_phase = 2*np.pi\n", + "phase_grid = np.linspace(min_phase,max_phase,n_phases)#between 0 - 2pi\n", + "#send params to phase angle routine\n", + "#case_3d.phase_angle(phase_grid=phase_grid,\n", + "# num_gangle=6, num_tangle=6,calculation='reflected')\n", + "\n", + "case_3d_clouds.phase_curve_geometry('reflected', phase_grid=phase_grid, num_gangle=6, num_tangle=6)\n", + "case_3d_clouds.inputs['atmosphere']['profile']\n", + "\n", + "case_3d_clouds.atmosphere_4d(gcm_out, shift = np.zeros(n_phases), zero_point='secondary_eclipse',\n", + " plot=True,verbose=False)\n", + "\n", + "# run clouds_4d (which must always run after atmosphere_4d)\n", + "case_3d_clouds.clouds_4d(ds_cld, iz_plot=34, plot=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Calculate cloudy phase curve**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "case_3d_clouds.gravity(radius=1,radius_unit=jdi.u.Unit('R_jup'),\n", + " mass=1, mass_unit=jdi.u.Unit('M_jup')) #any astropy units available\n", + "case_3d_clouds.star(opacity,5000,0,4.0, radius=1, radius_unit=jdi.u.Unit('R_sun'), semi_major=0.06, semi_major_unit=u.Unit('au'))\n", + "\n", + "allout_clouds = case_3d_clouds.phase_curve(opacity, n_cpu = 3,#jdi.cpu_count(),\n", + " full_output=False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "to_plot = 'fpfs_reflected' #or thermal\n", + "#collapse = [0.40 ]#micron\n", + "collapse='np.mean'\n", + "#phases, all_curves, all_ws, fig=jpi.phase_curve(allout, to_plot, collapse=collapse, R=R)\n", + "phases, all_curves, all_ws, fig=jpi.phase_curve(allout, to_plot, collapse=collapse, R=10,plot_width=1500, reorder_output=True)\n", + "phases_clouds, all_curves_clouds, all_ws_clouds, fig_clouds=jpi.phase_curve(allout_clouds, to_plot, collapse=collapse, R=10,plot_width=1500, reorder_output=True)\n", + "#jpi.show(fig)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compare cloudy and cloud free\n", + "\n", + "Below is a comparison between our cloudless and cloudy phase curves. The cloudy phase curve is much brighter. We also see a shift in the phase maximum due to our inhomogeneous cloud layer." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "plt.figure(figsize=(15, 7.5))\n", + "\n", + "plt.plot(phases,all_curves * 1e6, label=\"Cloudless\")\n", + "plt.plot(phases_clouds,all_curves_clouds * 1e6, label=\"Cloudy\")\n", + "\n", + "plt.xlabel('Orbital Phase')\n", + "plt.ylabel('F$_{P}$/F$_{S}$ (ppm)')\n", + "plt.title('Reflected Light Phase Curves')\n", + "plt.legend(loc='upper left')\n", + "plt.rcParams['font.size'] = 15\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Cloudy phase spectra\n", + "\n", + "The spectra for phase angles < 180 is flat due to clouds." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#same old same old\n", + "wno =[];fpfs=[];legend=[]\n", + "for iphase in allout_clouds.keys():\n", + " w,f = jdi.mean_regrid(allout_clouds[iphase]['wavenumber'],\n", + " allout_clouds[iphase]['fpfs_reflected'], R=100)\n", + " wno+=[w]\n", + " fpfs+=[f*1e6]\n", + " legend +=[str(int(iphase*180/np.pi))]\n", + "jpi.show(jpi.spectrum(wno, fpfs, plot_width=500,legend=legend,\n", + " palette=jpi.pals.viridis(n_phases)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To further analyze your output, you can use all of the same tools found in [Modeling a Phase Curve pt 1 (Robbins-Blanch et al. 2022) > Phase Curve Plotting Tools: Phase Snaps](https://natashabatalha.github.io/picaso/notebooks/9e_PhaseCurves.html#Phase-Curve-Plotting-Tools:-Phase-Snaps)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/docs/tutorials.rst b/docs/tutorials.rst index ca807676..7dfdebeb 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -51,8 +51,9 @@ Relevant Citatons: `Adams et al. 2022 Post-Processing Clouds for 3D runs Modeling a 3D Spectrum (Adams et al. 2022) - Modeling a Phase Curve pt 1 (Robbins-Blanch et al. 2022) - Modeling a Phase Curve pt 2 (Robbins-Blanch et al. 2022) + Modeling a Thermal Phase Curve pt 1 (Robbins-Blanch et al. 2022) + Modeling a Thermal Phase Curve pt 2 (Robbins-Blanch et al. 2022) + Modeling a Reflected Light Phase Curve (Hamill et al. 2024) 1D Climate Modeling ------------------- diff --git a/picaso/atmsetup.py b/picaso/atmsetup.py index d8336ba4..d557fa7f 100755 --- a/picaso/atmsetup.py +++ b/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: @@ -129,8 +134,11 @@ def get_profile_3d(self): if i == 'e-': electrons = True else: #don't raise exception, instead add user warning that a column has been automatically skipped - self.add_warnings("Ignoring %s in input file, not recognized molecule" % i) - warnings.warn("Ignoring %s in input file, not a recognized molecule" % i, UserWarning) + #make sure no warnin for new lat/lon keys + if 'lat' not in i: + if 'lon' not in i: + self.add_warnings("Ignoring %s in input file, not recognized molecule" % i) + warnings.warn("Ignoring %s in input file, not a recognized molecule" % i, UserWarning) first = False self.weights = weights diff --git a/picaso/disco.py b/picaso/disco.py index a8dbd762..46558beb 100755 --- a/picaso/disco.py +++ b/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 diff --git a/picaso/fluxes.py b/picaso/fluxes.py index 7c1e7f29..94f17420 100644 --- a/picaso/fluxes.py +++ b/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 diff --git a/picaso/justdoit.py b/picaso/justdoit.py index a3eb563b..e6a5ceec 100755 --- a/picaso/justdoit.py +++ b/picaso/justdoit.py @@ -21,6 +21,7 @@ from scipy.io import FortranFile + import os import pickle as pk import numpy as np @@ -442,6 +443,7 @@ def picaso(bundle,opacityclass, dimension = '1d',calculation='reflected', single_phase,multi_phase, frac_a,frac_b,frac_c,constant_back,constant_forward) xint_at_top += xint*gauss_wts[ig] + #if full output is requested add in xint at top for 3d plots if full_output: atm.xint_at_top = xint_at_top @@ -2671,12 +2673,17 @@ def atmosphere_4d(self, ds=None, shift=None, plot=True, iz_plot=0,verbose=True, if isinstance(shift, type(None)): shift = np.zeros(len(phases)) - if zero_point == 'night_transit': - shift = shift + 180 + if zero_point == 'night_transit': ## does not work for reflected case! + if 'reflected' in self.inputs['disco']['calculation']: + if verbose: print('Switching to zero point secondary_eclipse which is required for reflected light') + shift=shift + else: + if verbose: print('The zero_point input will be deprecated in the next PICASO version as it does not work for the reflectd light case. Instead things can be reordered in the phase_curve function in justplotit.phase_curve using reorder_output keyword') + shift = shift + 180 elif zero_point == 'secondary_eclipse': - shift=shift + shift=shift else: - raise Exception("Do not regocnize input zero point. Please specify: night_transit or secondary_eclipse") + raise Exception("Do not recognize input zero point. Please specify: night_transit or secondary_eclipse") self.inputs['shift'] = shift @@ -2690,7 +2697,6 @@ def atmosphere_4d(self, ds=None, shift=None, plot=True, iz_plot=0,verbose=True, #check for temperature and pressure if 'temperature' not in ds: raise Exception('Must include temperature as data component') - #check for pressure and change units if needed if 'pressure' not in ds.coords: raise Exception("Must include pressure in coords and units") @@ -2713,19 +2719,84 @@ def atmosphere_4d(self, ds=None, shift=None, plot=True, iz_plot=0,verbose=True, else : og_lat = ds.coords['lat'].values #degrees og_lon = ds.coords['lon'].values #degrees - + pres = ds.coords['pressure'].values #bars #adding this so I can add it explicitly to ds_New below + #store so we can rotate data_vars_og = {i:ds[i].values for i in ds.keys()} #run through phases and regrid each one + new_lat_totals = [] + new_lon_totals = [] + #new_lat_totals_og = [] + new_lon_totals_og = [] shifted_grids = {} + + # Add if calculation = reflected here: for i,iphase in enumerate(phases): new_lat = self.inputs['disco'][iphase]['latitude']*180/np.pi#to degrees - new_lon = self.inputs['disco'][iphase]['longitude']*180/np.pi#to degrees - total_shift = (iphase*180/np.pi + shift[i]) % 360 + new_lon_og = self.inputs['disco'][iphase]['longitude']*180/np.pi#to degrees + + + #Reflected case needs a step to ensure that the reflected crescent is at the correct point wrt the substellar point + #This statement only works for 10x10 cases! I am working on expanding this to all grids soon if possible + micro_shift = (abs(abs(new_lon_og[-1]) - abs(new_lon_og[-2])) - abs(abs(new_lon_og[0]) - abs(new_lon_og[1]))) / 2 #accounts for the difference in sizes between latxlon bins at phases!=0. + ng = self.inputs['disco'][iphase]['num_gangle'] # phase curve shift dependent on ngangles. Ng is used below to determine correct shift paramters + if ng == 6: + if 68 69 and new_lon_og[0] < 0 : #for first quarter of phases + new_lon_transfer = abs(new_lon_og[-1]) - abs(new_lon_og[0]) # take the difference between the first lon and the last lon at each phase + new_lon = new_lon_og - new_lon_transfer - micro_shift # The 'transfer' will then shift each phase to the opposite side of the dayside hemisphere. This is crucial for weighting ng and nt correctly for spectrum. + #add total shift statement + shift_back = -new_lon_transfer - micro_shift + elif new_lon_og[-1] > 69 and new_lon_og[0] > 0: # Second quarter of phases + new_lon_transfer = new_lon_og[-1] + new_lon_og[0] + new_lon = new_lon_og - new_lon_transfer - micro_shift + #print("new_lon 2nd Q", new_lon) + shift_back = -new_lon_transfer - micro_shift + elif new_lon_og[-1] < -69 and new_lon_og[0] < 0: # third quarter of phases + new_lon_transfer = new_lon_og[-1] + new_lon_og[0] + new_lon = new_lon_og - new_lon_transfer + micro_shift #new_lon_transfer here is negative, so we are adding + shift_back = -new_lon_transfer + micro_shift + elif new_lon_og[-1] < -69 and new_lon_og[0] > 0: #last quarter of phases + new_lon_transfer = abs(new_lon_og[-1]) - abs(new_lon_og[0]) # take the difference between the first lon and the last lon at each phase + new_lon = new_lon_og + new_lon_transfer + micro_shift# The 'transfer' will then shift each phase to the opposite side of the dayside hemisphere. This is crucial for weighting ng and nt correctly for spectrum. + #add total shift statement + shift_back = new_lon_transfer + micro_shift + if ng >= 10: + if 76 77 and new_lon_og[0] < 0 : #for first quarter of phases + new_lon_transfer = abs(new_lon_og[-1]) - abs(new_lon_og[0]) # take the difference between the first lon and the last lon at each phase + new_lon = new_lon_og - new_lon_transfer - micro_shift # The 'transfer' will then shift each phase to the opposite side of the dayside hemisphere. This is crucial for weighting ng and nt correctly for spectrum. + #add total shift statement + shift_back = -new_lon_transfer - micro_shift + elif new_lon_og[-1] > 77 and new_lon_og[0] > 0: # Second quarter of phases + new_lon_transfer = new_lon_og[-1] + new_lon_og[0] + new_lon = new_lon_og - new_lon_transfer - micro_shift + #print("new_lon 2nd Q", new_lon) + shift_back = -new_lon_transfer - micro_shift + elif new_lon_og[-1] < -77 and new_lon_og[0] < 0: # third quarter of phases + new_lon_transfer = new_lon_og[-1] + new_lon_og[0] + new_lon = new_lon_og - new_lon_transfer + micro_shift #new_lon_transfer here is negative, so we are adding + shift_back = -new_lon_transfer + micro_shift + elif new_lon_og[-1] < -77 and new_lon_og[0] > 0: #last quarter of phases + new_lon_transfer = abs(new_lon_og[-1]) - abs(new_lon_og[0]) # take the difference between the first lon and the last lon at each phase + new_lon = new_lon_og + new_lon_transfer + micro_shift# The 'transfer' will then shift each phase to the opposite side of the dayside hemisphere. This is crucial for weighting ng and nt correctly for spectrum. + #add total shift statement + shift_back = new_lon_transfer + micro_shift + + #append new lons and lats, used to create array below this for loop + new_lat_totals.append(new_lat) + new_lon_totals.append(new_lon) + new_lon_totals_og.append(new_lon_og) + self.inputs['disco'][iphase]['longitude'] = new_lon * np.pi/180 # changing lon around requires us to re-define self.inputs as well (needed for disco geom and also for clouds_4d) + total_shift = (iphase*180/np.pi + (shift[i] + shift_back)) % 360 + #total_shift = (iphase*180/np.pi + shift[i]) % 360 change_zero_pt = og_lon + total_shift change_zero_pt[change_zero_pt>360]=change_zero_pt[change_zero_pt>360]%360 #such that always between -180 and 180 change_zero_pt[change_zero_pt>180]=change_zero_pt[change_zero_pt>180]%180-180 #such that always between -180 and 180 - #ds.coords['lon'].values = change_zero_pt split = np.argmin(abs(change_zero_pt + 180)) #find point where we should shift the grid for idata in data_vars_og.keys(): swap1 = data_vars_og[idata][0:split,:,:] @@ -2733,14 +2804,45 @@ def atmosphere_4d(self, ds=None, shift=None, plot=True, iz_plot=0,verbose=True, data = np.concatenate((swap2,swap1)) ds[idata].values = data shifted_grids[iphase] = regrid_xarray(ds, latitude=new_lat, longitude=new_lon) - new_phase_grid=xr.concat(list(shifted_grids.values()), pd.Index(list(shifted_grids.keys()), name='phase')) - + + # we need arrays that are len(phase) x len(lon regrid) as array, not list. + # These are used to create 'lon2d' and 'lat2d', which are needed for reflected case. + new_lat_totals_array = np.array(new_lat_totals) + new_lon_totals_array = np.array(new_lon_totals) + new_lon_totals_og_array = np.array(new_lon_totals_og) + + # This creates 'phase' as a coord + stacked_phase_grid=xr.concat(list(shifted_grids.values()), pd.Index(list(shifted_grids.keys()), name='phase'), join='override') ## join=override gets rid of errant lon values + + # Here we are manually creating a new xarray from scratch that has 'lon2d', 'lat2d', which have 'phase' as their 2nd dimension (neeeded for reflected case) + # This is a temporary xarray that will be used to merge data variables (created above) with our new 2d coordinates. + # We do it this way because xarray does not like when you add dimensions to existing coordinate system. This seems to be the only work around. + ds_New = xr.Dataset( + data_vars=dict( + ), + coords=dict( + lon2d=(["phase","lon"], new_lon_totals_array,{'units': 'degrees'}), #required. Errors when named lon + lat2d=(["phase","lat"], new_lat_totals_array,{'units': 'degrees'}), #required + lon2d_clouds=(["phase","lon"], new_lon_totals_og_array,{'units': 'degrees'}), #This is the original coord system. We need to conserve this for clouds_4d + lat2d_clouds=(["phase","lat"], new_lat_totals_array,{'units': 'degrees'}), # lon2d_clouds and lat2d_clouds will be used for shift in clouds_4d. This is the only purpose of these two coords. + pressure=(["pressure"], pres,{'units': 'bar'}) #required + ), + attrs=dict(description="coords with vectors"), + ) + new_phase_grid = ds_New + + # Lets use merge with compat=override (use data_vars from 1st dataset) and join=right + # This creates an xarray with all of the variables from stacked_phase_grid (i.e., temperature and chemicals). + # This also creates an xarray with coords named 'lon2d' and 'lat2d' (as well as 'lon' and 'lat'). 'lon2d' and 'lat2d' have 'phase' as their second dimension, which is needed when we use reflected case. + new_phase_grid = xr.merge([stacked_phase_grid, new_phase_grid], compat='override', join='right') + if plot: - new_phase_grid['temperature'].isel(pressure=iz_plot).plot(x='lon', y ='lat', col='phase',col_wrap=4) + new_phase_grid['temperature'].isel(pressure=iz_plot).plot(x='lon2d', y ='lat2d', col='phase',col_wrap=4) + #changed lon, lat to lon2d, lat2d self.inputs['atmosphere']['profile'] = new_phase_grid - def clouds_4d(self, ds=None, plot=True, iz_plot=0,iw_plot=0,verbose=True): + def clouds_4d(self, ds=None, plot=True, iz_plot=0,iw_plot=0,verbose=True, calculation='reflected'): """ Regrids xarray @@ -2781,7 +2883,6 @@ def clouds_4d(self, ds=None, plot=True, iz_plot=0,iw_plot=0,verbose=True): if 'g0' not in ds: raise Exception('Must include g0 as data component') if 'w0' not in ds: raise Exception('Must include w0 as data component') - #check for pressure and change units if needed if 'pressure' not in ds.coords: raise Exception("Must include pressure in coords and units") @@ -2816,33 +2917,147 @@ def clouds_4d(self, ds=None, plot=True, iz_plot=0,iw_plot=0,verbose=True): else : og_lat = ds.coords['lat'].values #degrees og_lon = ds.coords['lon'].values #degrees + pres = ds.coords['pressure'].values #bars #adding this so I can add it explicitly to ds_New below + if 'reflected' in calculation: + #store so we can rotate + data_vars_og = {i:ds[i].values for i in ds.keys()} + #run through phases and regrid each one + new_lat_totals = [] + new_lon_totals = [] + shifted_grids = {} + for i,iphase in enumerate(phases): + new_lat = np.array(self.inputs['atmosphere']['profile']['lat2d_clouds'][i,:])#*180/np.pi + new_lon_og = np.array(self.inputs['atmosphere']['profile']['lon2d_clouds'][i,:])#*180/np.pi + micro_shift = (abs(abs(new_lon_og[-1]) - abs(new_lon_og[-2])) - abs(abs(new_lon_og[0]) - abs(new_lon_og[1]))) / 2 #accounts for the difference in sizes between latxlon bins at phases!=0. + ng = self.inputs['disco'][iphase]['num_gangle'] # phase curve shift dependent on ngangles. Ng is used below to determine correct shift paramters + if ng == 6: + if 68 69 and new_lon_og[0] < 0 : #for first quarter of phases + new_lon_transfer = abs(new_lon_og[-1]) - abs(new_lon_og[0]) # take the difference between the first lon and the last lon at each phase + new_lon = new_lon_og - new_lon_transfer - micro_shift # The 'transfer' will then shift each phase to the opposite side of the dayside hemisphere. This is crucial for weighting ng and nt correctly for spectrum. + #add total shift statement + shift_back = -new_lon_transfer - micro_shift + elif new_lon_og[-1] > 69 and new_lon_og[0] > 0: # Second quarter of phases + new_lon_transfer = new_lon_og[-1] + new_lon_og[0] + new_lon = new_lon_og - new_lon_transfer - micro_shift + #print("new_lon 2nd Q", new_lon) + shift_back = -new_lon_transfer - micro_shift + elif new_lon_og[-1] < -69 and new_lon_og[0] < 0: # third quarter of phases + new_lon_transfer = new_lon_og[-1] + new_lon_og[0] + new_lon = new_lon_og - new_lon_transfer + micro_shift #new_lon_transfer here is negative, so we are adding + shift_back = -new_lon_transfer + micro_shift # - 180 + elif new_lon_og[-1] < -69 and new_lon_og[0] > 0: #last quarter of phases + new_lon_transfer = abs(new_lon_og[-1]) - abs(new_lon_og[0]) # take the difference between the first lon and the last lon at each phase + new_lon = new_lon_og + new_lon_transfer + micro_shift # The 'transfer' will then shift each phase to the opposite side of the dayside hemisphere. This is crucial for weighting ng and nt correctly for spectrum. + #add total shift statement + shift_back = new_lon_transfer + micro_shift # - 180 + if ng >= 10: + if 76 77 and new_lon_og[0] < 0 : #for first quarter of phases + new_lon_transfer = abs(new_lon_og[-1]) - abs(new_lon_og[0]) # take the difference between the first lon and the last lon at each phase + new_lon = new_lon_og - new_lon_transfer - micro_shift # The 'transfer' will then shift each phase to the opposite side of the dayside hemisphere. This is crucial for weighting ng and nt correctly for spectrum. + #add total shift statement + shift_back = -new_lon_transfer - micro_shift + elif new_lon_og[-1] > 77 and new_lon_og[0] > 0: # Second quarter of phases + new_lon_transfer = new_lon_og[-1] + new_lon_og[0] + new_lon = new_lon_og - new_lon_transfer - micro_shift + #print("new_lon 2nd Q", new_lon) + shift_back = -new_lon_transfer - micro_shift + elif new_lon_og[-1] < -77 and new_lon_og[0] < 0: # third quarter of phases + new_lon_transfer = new_lon_og[-1] + new_lon_og[0] + new_lon = new_lon_og - new_lon_transfer + micro_shift #new_lon_transfer here is negative, so we are adding + shift_back = -new_lon_transfer + micro_shift - 180 + elif new_lon_og[-1] < -77 and new_lon_og[0] > 0: #last quarter of phases + new_lon_transfer = abs(new_lon_og[-1]) - abs(new_lon_og[0]) # take the difference between the first lon and the last lon at each phase + new_lon = new_lon_og + new_lon_transfer + micro_shift # The 'transfer' will then shift each phase to the opposite side of the dayside hemisphere. This is crucial for weighting ng and nt correctly for spectrum. + #add total shift statement + shift_back = new_lon_transfer + micro_shift - 180 + + new_lat_totals.append(new_lat) + new_lon_totals.append(new_lon) + #total_shift = (iphase*180/np.pi + (shift[i] - shift_back)) % 360 + # total_shift = (iphase*180/np.pi + (shift[i] + shift_back)) % 360 + total_shift = (iphase*180/np.pi + shift[i]) % 360 + change_zero_pt = og_lon + total_shift + shift_back + change_zero_pt[change_zero_pt>360]=change_zero_pt[change_zero_pt>360]%360 #such that always between -180 and 180 + change_zero_pt[change_zero_pt>180]=change_zero_pt[change_zero_pt>180]%180-180 #such that always between -180 and 180 + split = np.argmin(abs(change_zero_pt + 180)) #find point where we should shift the grid + for idata in data_vars_og.keys(): + swap1 = data_vars_og[idata][0:split,:,:] + swap2 = data_vars_og[idata][split:,:,:] + data = np.concatenate((swap2,swap1)) + ds[idata].values = data + shifted_grids[iphase] = regrid_xarray(ds, latitude=new_lat, longitude=new_lon) + + ## we need a lon_total that is len(phase) x len(lon regrid) as ARRAY, not list + new_lat_totals_array = np.array(new_lat_totals) + new_lon_totals_array = np.array(new_lon_totals) + + # creates phase as a coord + stacked_phase_grid=xr.concat(list(shifted_grids.values()), pd.Index(list(shifted_grids.keys()), name='phase'), join='override') ## join=override gets rid of errant lon values + + # put data into a dataset + ds_New = xr.Dataset( + data_vars=dict( + ), + coords=dict( + lon2d=(["phase","lon"], new_lon_totals_array,{'units': 'degrees'}), #required. Errors when named lon + lat2d=(["phase","lat"], new_lat_totals_array,{'units': 'degrees'}), #required + pressure=(["pressure"], pres,{'units': 'bar'})#required* + ), + attrs=dict(description="coords with vectors"), + ) + new_phase_grid = ds_New - #store so we can rotate - data_vars_og = {i:ds[i].values for i in ds.keys()} - #run through phases and regrid each one - shifted_grids = {} - for i,iphase in enumerate(phases): - new_lat = self.inputs['disco'][iphase]['latitude']*180/np.pi#to degrees - new_lon = self.inputs['disco'][iphase]['longitude']*180/np.pi#to degrees - total_shift = (iphase*180/np.pi + shift[i]) % 360 - change_zero_pt = og_lon + total_shift - change_zero_pt[change_zero_pt>360]=change_zero_pt[change_zero_pt>360]%360 #such that always between -180 and 180 - change_zero_pt[change_zero_pt>180]=change_zero_pt[change_zero_pt>180]%180-180 #such that always between -180 and 180 - #ds.coords['lon'].values = change_zero_pt - split = np.argmin(abs(change_zero_pt + 180)) #find point where we should shift the grid - for idata in data_vars_og.keys(): - swap1 = data_vars_og[idata][0:split,:,:] - swap2 = data_vars_og[idata][split:,:,:] - data = np.concatenate((swap2,swap1)) - ds[idata].values = data - shifted_grids[iphase] = regrid_xarray(ds, latitude=new_lat, longitude=new_lon) - new_phase_grid=xr.concat(list(shifted_grids.values()), pd.Index(list(shifted_grids.keys()), name='phase')) + # Now we need to add stacked_phase_grid Data Vars to new_phase_grid, and also add Phase to coords + + # Lets use merge with compat=override (use data_vars from 1st dataset) + # This adds a new, 2D coord named 'lon2d' (not 'lon') and 'lat2d' (not 'lon). Lon2d needs to be specified for phase_curve + new_phase_grid = xr.merge([stacked_phase_grid, new_phase_grid], compat='override', join='right') - if plot: - new_phase_grid['opd'].isel(pressure=iz_plot,wno=iw_plot).plot(x='lon', y ='lat', col='phase',col_wrap=4) - - self.inputs['clouds']['profile'] = new_phase_grid - self.inputs['clouds']['wavenumber'] = ds.coords['wno'].values + #print(" Cloud Phase Grid XArray", new_phase_grid) + + if plot: + #new_phase_grid['opd'].isel(pressure=iz_plot,wno=iw_plot).plot(x='lon2d', y ='lat2d', col='phase',col_wrap=4) + new_phase_grid['opd'].isel(pressure=iz_plot,wno=iw_plot).plot(x='lon2d', y ='lat2d', col='phase',col_wrap=4) + + self.inputs['clouds']['profile'] = new_phase_grid + self.inputs['clouds']['wavenumber'] = ds.coords['wno'].values + + elif 'thermal' in calculation: # copy-paste of original clouds_4d + #store so we can rotate + data_vars_og = {i:ds[i].values for i in ds.keys()} + #run through phases and regrid each one + shifted_grids = {} + for i,iphase in enumerate(phases): + new_lat = self.inputs['disco'][iphase]['latitude']*180/np.pi#to degrees + new_lon = self.inputs['disco'][iphase]['longitude']*180/np.pi#to degrees + total_shift = (iphase*180/np.pi + shift[i]) % 360 + change_zero_pt = og_lon + total_shift + change_zero_pt[change_zero_pt>360]=change_zero_pt[change_zero_pt>360]%360 #such that always between -180 and 180 + change_zero_pt[change_zero_pt>180]=change_zero_pt[change_zero_pt>180]%180-180 #such that always between -180 and 180 + #ds.coords['lon'].values = change_zero_pt + split = np.argmin(abs(change_zero_pt + 180)) #find point where we should shift the grid + for idata in data_vars_og.keys(): + swap1 = data_vars_og[idata][0:split,:,:] + swap2 = data_vars_og[idata][split:,:,:] + data = np.concatenate((swap2,swap1)) + ds[idata].values = data + shifted_grids[iphase] = regrid_xarray(ds, latitude=new_lat, longitude=new_lon) + new_phase_grid=xr.concat(list(shifted_grids.values()), pd.Index(list(shifted_grids.keys()), name='phase')) + + if plot: + new_phase_grid['opd'].isel(pressure=iz_plot,wno=iw_plot).plot(x='lon', y ='lat', col='phase',col_wrap=4) + + self.inputs['clouds']['profile'] = new_phase_grid + self.inputs['clouds']['wavenumber'] = ds.coords['wno'].values + + else: + raise Exception("Must include 'reflected' or 'thermal' in calculation") def surface_reflect(self, albedo, wavenumber, old_wavenumber = None): """ @@ -3148,6 +3363,7 @@ def run_virga(ilon,ilat): results = Parallel(n_jobs=n_cpu)(delayed(run_virga)(ilon,ilat) for ilon in range(ng) for ilat in range(nt)) wno_grid = 1e4/results[0]['wave'] + wno_grid_sorted = sorted(wno_grid) nwno = len(wno_grid) pres = results[0]['pressure'] @@ -3173,7 +3389,8 @@ def run_virga(ilon,ilat): lon=(["lon"], lon,{'units': 'degrees'}),#required lat=(["lat"], lat,{'units': 'degrees'}),#required pressure=(["pressure"], pres,{'units': 'bar'}),#required - wno=(["wno"], wno_grid,{'units': 'cm^(-1)'})#required for clouds + # wno=(["wno"], wno_grid,{'units': 'cm^(-1)'})#required for clouds + wno=(["wno"], wno_grid_sorted,{'units': 'cm^(-1)'})#required for clouds ), attrs=dict(description="coords with vectors"), ) @@ -3417,7 +3634,7 @@ def approx(self,single_phase='TTHG_ray',multi_phase='N=2',delta_eddington=True, def phase_curve(self, opacityclass, full_output=False, - plot_opacity= False,n_cpu =1 ): + plot_opacity= False,n_cpu =1,verbose=True ): """ Run phase curve Parameters @@ -3433,12 +3650,16 @@ def phase_curve(self, opacityclass, full_output=False, phases = self.inputs['phase_angle'] calculation = self.inputs['disco']['calculation'] all_geom = self.inputs['disco'] + #print("all_geom", all_geom) all_profiles = self.inputs['atmosphere']['profile'] all_cld_profiles = self.inputs['clouds']['profile'] def run_phases(iphase): self.inputs['phase_angle'] = iphase[1] self.inputs['atmosphere']['profile'] = all_profiles.isel(phase=iphase[0]) + + if verbose: print("Currently computing Phase", iphase) + self.inputs['disco'] = all_geom[iphase[1]] if not isinstance(all_cld_profiles, type(None)): self.inputs['clouds']['profile'] = all_cld_profiles.isel(phase=iphase[0]) diff --git a/picaso/justplotit.py b/picaso/justplotit.py index 2eb2142d..f9d1283f 100755 --- a/picaso/justplotit.py +++ b/picaso/justplotit.py @@ -1559,7 +1559,10 @@ def phase_snaps(allout, x = 'longitude', y = 'pressure', z='temperature',palette fig.tight_layout() return fig -def phase_curve(allout, to_plot, collapse=None, R=100, palette=pals.Spectral11,verbose=True, **kwargs): + +def phase_curve(allout, to_plot, collapse=None, R=100, + palette=pals.Spectral11,verbose=True, + reorder_output=False, **kwargs): """ Plots phase curves @@ -1581,7 +1584,11 @@ def phase_curve(allout, to_plot, collapse=None, R=100, palette=pals.Spectral11,v Print out low level warnings kwargs : dict Bokeh plotting kwargs for bokeh.Figure + reorder_output : bool + Returns sorted outputs, for better phase curve plotting. + re-orders phases such that brightest value is at phase=0 """ + kwargs['height'] = kwargs.get('plot_height',kwargs.get('height',400)) kwargs['width'] = kwargs.get('plot_width', kwargs.get('width',600)) if 'plot_width' in kwargs.keys() : kwargs.pop('plot_width') @@ -1632,10 +1639,46 @@ def phase_curve(allout, to_plot, collapse=None, R=100, palette=pals.Spectral11,v legend = Legend(items=legend_it, location=(0, -20)) legend.click_policy="mute" fig.add_layout(legend, 'left') - + fig.xgrid.grid_line_alpha=0 fig.ygrid.grid_line_alpha=0 plot_format(fig) + + #if verbose: print("phases", phases) + #if verbose: print("all_curves",all_curves) + + #re-order phases such that brightest value is at phase=0 (only way to do this for reflected case) + front_half_phases = phases[:len(phases)//2] + #print("front_half_phases", front_half_phases) + back_half_phases = phases[len(phases)//2:] - (2*np.pi) + #print("back_half_phases", back_half_phases) + reorder_phases = np.concatenate((back_half_phases, front_half_phases)) + + #re-order all_curves such that brightest value is at phase=0 (only way to do this for reflected case) + front_half_all_curves = all_curves[:len(all_curves)//2] + back_half_all_curves = all_curves[len(all_curves)//2:] + reorder_all_curves = np.concatenate((back_half_all_curves, front_half_all_curves)) + #print("Reorder all curves",reorder_all_curves) + + if to_plot == "fpfs_reflected" or to_plot == "albedo": + fig2 = figure(**kwargs) + + for i in range(len(collapse)): + f2 = fig2.line(reorder_phases,reorder_all_curves[:,i],line_width=3,color=palette[i], + ) + + legend2 = Legend(items=legend_it, location=(0, -20)) + legend2.click_policy="mute" + fig2.add_layout(legend, 'left') + + fig2.xgrid.grid_line_alpha=0 + fig2.ygrid.grid_line_alpha=0 + plot_format(fig2) + #show(fig2) + + if reorder_output: + return reorder_phases, reorder_all_curves, all_ws, fig + return phases, all_curves, all_ws, fig def thermal_contribution(full_output, tau_max=1.0,R=100, **kwargs):