plot_template_Slab89.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mgimg\n",
    "import matplotlib.colors as colors\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "import emcee\n",
    "import testkit\n",
    "import corner\n",
    "import pickle as pickle\n",
    "from IPython.display import display\n",
    "%matplotlib inline\n",
    "import forwardmodel\n",
    "import ciamod\n",
    "import TPmod\n",
    "import cloud\n",
    "import band\n",
    "import brewtools\n",
    "from astropy.convolution import convolve, convolve_fft\n",
    "from astropy.convolution import Gaussian1DKernel\n",
    "from scipy.io.idl import readsav\n",
    "from scipy import interpolate\n",
    "from scipy.interpolate import interp1d\n",
    "from scipy.interpolate import InterpolatedUnivariateSpline\n",
    "from bensconv import spex_non_uniform\n",
    "from bensconv import conv_uniform_R\n",
    "from bensconv import conv_uniform_FWHM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First step is to load the results file and the run arguments.\n",
    "\n",
    "We open the files and select the last couple of thousand iterations of the chain, and flatted this into a simple array of state vectors called \"samples\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: '<yourrunname>_snapshot.pic'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-647bf553196a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m# OK finish? 1 for yes, 0 for no.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mfin\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mflatendchain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflatendprobs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbrewtools\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_endchain\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrunname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfin\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0mtheta_max_end\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mflatendchain\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mflatendprobs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0mmax_end_like\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mamax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mflatendprobs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Documents/BDNYC/Retrievals/brewster/brewtools.py\u001b[0m in \u001b[0;36mget_endchain\u001b[0;34m(runname, fin, results_path)\u001b[0m\n\u001b[1;32m     24\u001b[0m     \u001b[0;32melif\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfin\u001b[0m \u001b[0;34m==\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     25\u001b[0m         \u001b[0mpic\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mresults_path\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mrunname\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m\"_snapshot.pic\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 26\u001b[0;31m         \u001b[0mchain\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mprobs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpickle_load\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpic\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     27\u001b[0m         \u001b[0mnwalkers\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mchain\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     28\u001b[0m         \u001b[0mntot\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mchain\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Documents/BDNYC/Retrievals/brewster/brewtools.py\u001b[0m in \u001b[0;36mpickle_load\u001b[0;34m(file_path)\u001b[0m\n\u001b[1;32m    227\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    228\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mpickle_load\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfile_path\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 229\u001b[0;31m     \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfile_path\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"rb\"\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    230\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mpickle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mMacOSFile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: '<yourrunname>_snapshot.pic'"
     ]
    }
   ],
   "source": [
    "#path = \"/you/can/add/a/path/to/your/results/files\"\n",
    "runname = \"<yourrunname>\"\n",
    "# OK finish? 1 for yes, 0 for no.\n",
    "fin = 0\n",
    "flatendchain, flatendprobs,ndim = brewtools.get_endchain(runname,fin)\n",
    "theta_max_end = flatendchain[np.argmax(flatendprobs)]\n",
    "max_end_like = np.amax(flatendprobs)\n",
    "samples = flatendchain\n",
    "\n",
    "argfile =runname+\"_runargs.pic\"\n",
    "\n",
    "runargs = brewtools.pickle_load(argfile)\n",
    "# If you're opening on a Linux box, use the code below\n",
    "#with open(argfile, 'rb') as input:\n",
    "#    runargs = pickle.load(input) \n",
    "\n",
    "gases_myP,chemeq,dist, cloudtype,do_clouds,gasnum,cloudnum,inlinetemps,coarsePress,press,inwavenum,linelist,cia,ciatemps,use_disort,fwhm,obspec,proftype,do_fudge, prof,do_bff,bff_raw,ceTgrid,metscale,coscale = runargs\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we print the max likelihood value for the state vector, just to have a look at it\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(theta_max_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And print the BIC for this run. You might want to take a note of this.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BIC = (-2.* max_end_like) + (ndim*np.log(obspec.shape[1]))\n",
    "print(\"BIC = \"+str(BIC))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next we want to plot the profile. \n",
    "\n",
    "We start by calculating the profiles from the temperature parameters in the samples array.\n",
    "\n",
    "This case assumes a type 2 profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tsamples = samples[:,ndim-5:]\n",
    "nsamps = Tsamples.shape[0]\n",
    "Tprofs = np.empty([64,Tsamples.shape[0]])\n",
    "for i in range(0,nsamps):\n",
    "    Tprofs[:,i] = TPmod.set_prof(2,coarsePress,press,Tsamples[i,:])\n",
    "    \n",
    "Tlays = np.empty([64,5])\n",
    "for i in range(0,64):\n",
    "    junk = Tprofs[i,:]\n",
    "    junk2 = np.percentile(junk, [2.4,16, 50, 84,97.6],axis=0)\n",
    "    junk3 = np.array(junk2)\n",
    "    Tlays[i,:] = junk3[:]\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to sort out the cloud parameters. this next cell is just about unpacking the cloud parameters from the state vector "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the cloud stuff\n",
    "cloudtemp = samples[:,ndim-10:ndim-5]\n",
    "cloudsamples = np.empty_like(cloudtemp)\n",
    "print(cloudtemp.shape)\n",
    "cloudsamples[:,0] = cloudtemp[:,0]\n",
    "cloudsamples[:,1] = cloudtemp[:,1] + cloudtemp[:,2]\n",
    "cloudsamples[:,2] = cloudtemp[:,2]\n",
    "cloudsamples[:,3] = cloudtemp[:,3]\n",
    "cloudsamples[:,4] = cloudtemp[:,4]\n",
    "\n",
    "cloud_Base = np.percentile(cloudsamples[:, 1], [16, 50, 84])\n",
    "cloud_height = np.percentile(cloudsamples[:, 2], [16, 50, 84])\n",
    "cloud_Tau = np.percentile(cloudsamples[:, 0], [16, 50, 84])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot it..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.rc('font',family='Times New Roman')\n",
    "fig=plt.figure(dpi=320)\n",
    "plt.axis([0., 4000.,3.0,-5.0])\n",
    "\n",
    "logP = np.log10(press)\n",
    "s = 'cloud_tau = '+str(round(cloud_Tau[0],2))+', '+str(round(cloud_Tau[1],2))+', '+str(round(cloud_Tau[2],2))\n",
    "\n",
    "plt.fill_between([500.,3500.0],cloud_Base[1],cloud_Base[1] - cloud_height[1],facecolor='blue',alpha=0.3)\n",
    "plt.text(500.,2.5,s)\n",
    "plt.fill_between([0.,500.0],cloud_Base[2],cloud_Base[2] - cloud_height[0],facecolor='grey',alpha=0.5)\n",
    "plt.fill_between([3500.,4000.0],cloud_Base[0],cloud_Base[0] - cloud_height[2],facecolor='grey',alpha=0.5)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "d1, = plt.plot(Tlays[:,2],logP,'k-',label=runname)\n",
    "plt.fill_betweenx(logP,Tlays[:,1], Tlays[:,3], facecolor='red', alpha=0.3)\n",
    "plt.fill_betweenx(logP,Tlays[:,0], Tlays[:,4], facecolor='red', alpha=0.1)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "enst = 10.0**4/(6.26 - 0.35*logP-0.70*0.0)\n",
    "fost = 10.0**4/(5.89 - 0.37*logP-0.73*0.0)\n",
    "iron = 10.0**4/(5.44 - 0.48*logP-0.48*0.0)\n",
    "cr =  10.0**4/(6.528 - 0.491*logP-0.491*0.0)\n",
    "al2o3 = 10.0**4 / (5.0139 - 0.21794*(logP) + 2.2636E-03*(logP)**2.0 - 0.580*0.0)\n",
    "c1, = plt.plot(enst,logP,'--',color='blue',linewidth=1.5, label='MgSiO$_3$')\n",
    "c2, = plt.plot(fost,logP,'--',color='pink',linewidth=1.5,label='Mg$_2$SiO$_4$')\n",
    "c3, = plt.plot(iron,logP,'--',color='orange',linewidth=1.5, label='Fe')\n",
    "c4, = plt.plot(cr,logP,'--',color='purple',linewidth=1.5, label='Cr')\n",
    "c5, = plt.plot(al2o3,logP,'--',color='red',linewidth=1.5, label='Al$_2$O$_3$')\n",
    "\n",
    "plt.legend(handles=[d1,c1,c2,c3,c5])\n",
    "plt.ylabel(r'log(P / bar)')\n",
    "plt.xlabel('T / K')\n",
    "\n",
    "plt.savefig(runname+\"_profile.png\",format='png', dpi=320)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Save the median profile\n",
    "np.savetxt(runname+'_profile.dat', np.c_[logP, Tlays[:,2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next we want to plot the gas abundances and cloud parameters in corner plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# You'll need to edit this cell to make it work for the gases you've used\n",
    "# e.g. change the cut you take out of the samples array to get gases and gravity\n",
    "gassamples = samples[:,0:10]\n",
    "# set up boundaries for histograms\n",
    "rh2o = (-12.0,.0)\n",
    "rco = (-12.0,0.)\n",
    "rco2 = (-12.0, 0.)\n",
    "rch4 = (-12,0.)\n",
    "rvo= (-12,0.)\n",
    "rtio = (-12,-0)\n",
    "#rcah = (-12,-0)\n",
    "rcrh = (-12,0)\n",
    "rfeh = (-12.0,-0.0)\n",
    "rnak = (-12.,0)\n",
    "rlogg = (3.5,5.5)\n",
    "# now make an array of the bounds to give to corner plot\n",
    "bnds = [rh2o,rco,rco2,rch4,rtio,rvo,rcrh,rfeh,rnak,rlogg]\n",
    "\n",
    "fig = corner.corner(gassamples,scale_hist=False, range = bnds,plot_datapoints =False,labels=[\"H2O\",\"CO\",\"CO2\",\"CH4\",\"TiO\",\"VO\",\"CrH\",\"FeH\",\"Na+K\",\"logg\"],quantiles=[0.16, 0.5, 0.84],show_titles=True, title_kwargs={\"fontsize\": 20},label_kwargs={\"fontsize\": 20})\n",
    "#plt.savefig(runname+\"_gascorner.png\",format='png', dpi=320)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This saves out the gas mixing ratios which might be needed later\n",
    "gasVMR = np.ones([gasnum.size-1,3])\n",
    "for i in range (0,gasnum.size-1):\n",
    "    gasVMR[i,:] = np.percentile(samples[:,i], [16, 50, 84])\n",
    "    \n",
    "np.savetxt(runname+'_VMRs.dat', np.c_[gasVMR[:,0],gasVMR[:,1],gasVMR[:,2]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This plots a corner plot for the cloud parameters\n",
    "rtau1 = (1,8)\n",
    "rp01 = (0.,2.0)\n",
    "rdp1 = (0,3)\n",
    "rw0 = (0,1)\n",
    "rpow = (-3,+3)\n",
    "\n",
    "bnds = [rtau1,rp01,rdp1,rw0,rpow]\n",
    "fig = corner.corner(cloudsamples,range =bnds,scale_hist=False,plot_datapoints =False,labels=[\"tau1\",\"$log(P_{base}  / bar)$\",\"dP\",\"w0\",\"pow\"],quantiles=[0.16, 0.5, 0.84],show_titles=True, title_kwargs={\"fontsize\": 20,\"loc\":'left'},label_kwargs={\"fontsize\": 20})\n",
    "plt.savefig(runname+\"_cloudcorner.png\",format='png', dpi=320)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting the model spectra\n",
    "To plot the model spectra we need to rerun the model, and rerun it for a bunch of random draws from the posterior. This will allow us to plot either a spaghetti plot, or a median + interval spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get diagnostics along with the spectrum\n",
    "gnostics = 0\n",
    "# Now run the model again to get your model spectrum and process to make it look like the data\n",
    "shiftspec, photspec, tauspec,cfunc = testkit.modelspec(theta_max_end,runargs,gnostics)\n",
    "topspec = brewtools.proc_spec(shiftspec,theta_max_end,fwhm,chemeq,gasnum,obspec) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now grab 500 random draws from the posterior\n",
    "pltspec = np.zeros((500,obspec[0,:].size))\n",
    "samp= np.empty(ndim)\n",
    "samples = flatendchain\n",
    "sid = np.zeros(500)\n",
    "for i in range (0,500):\n",
    "    sid[i]= np.random.randint(0,high = len(samples))\n",
    "    samp = samples[int(sid[i]),:]\n",
    "    shiftspec, photspec, tauspec,cfunc = testkit.modelspec(samp,runargs,gnostics)\n",
    "    pltspec[i,:] = brewtools.proc_spec(shiftspec,samp,fwhm,chemeq,gasnum,obspec) \n",
    "\n",
    "# get the intervals for the distribution of model spectra\n",
    "specdist = np.empty([obspec[0].size,5])\n",
    "for i in range(0,obspec[0].size):\n",
    "    junk = pltspec[:,i]\n",
    "    junk2 = np.percentile(junk, [2.4,16, 50, 84,97.6],axis=0)\n",
    "    junk3 = np.array(junk2)\n",
    "    specdist[i,:] = junk3[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the spectra\n",
    "plt.rc('font',family='Times New Roman')\n",
    "fig=plt.figure(dpi=320)\n",
    "plt.axis([1.0,2.5,0,1e-14])\n",
    "# you can swap in the lines below for a log plot\n",
    "#plt.axis([0.6,15.0, 1.e-17,150e-15])\n",
    "\n",
    "#ax = fig.add_subplot(1,1,1)\n",
    "\n",
    "#ax.set_yscale('log')\n",
    "#ax.set_xscale('log')\n",
    "\n",
    "# add in below lines for spaghetti plot\n",
    "#for i in range(0,500):\n",
    "#    plt.plot(obspec[0],pltspec[i,:],color='r',alpha=0.01)\n",
    "\n",
    "d1, = plt.plot(obspec[0,:],obspec[1,:],'k-',label = \"data\")\n",
    "t1, = plt.plot(obspec[0,:],topspec,'g-',linewidth=1, label = runname+\" max likelihood\")\n",
    "\n",
    "r1, = plt.plot(obspec[0],specdist[:,2],'y-',linewidth=0.5, label = \"median\")\n",
    "#plt.fill_between(obspec[0],specdist[:,0],specdist[:,4],facecolor='red',alpha=0.2)\n",
    "plt.fill_between(obspec[0],specdist[:,1],specdist[:,3],facecolor='red',alpha=0.5)\n",
    "\n",
    "#plt.fill_between(obspec[0,:],obspec[1,:]-obspec[2,:],obspec[1,:]+obspec[2,:],facecolor='red',alpha=0.2)\n",
    "\n",
    "\n",
    "plt.legend(handles=[d1,t1])\n",
    "\n",
    "\n",
    "plt.ylabel(r'$ F_{\\lambda}$ / $Wm^{-2} \\mu m^{-1}$')\n",
    "plt.xlabel('Wavelength / $\\mu m$')\n",
    "#plt.savefig(runname+\"_SPAG_SPEC.png\",format='png', dpi=320)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Contribution function\n",
    "Finally, the contribution function. This is a really useful plot..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This one plots a contribution function \n",
    "\n",
    "# get diagnostics along with the spectrum\n",
    "gnostics = 1\n",
    "shiftspec, clphotspec, ophotspec,cfunc = testkit.modelspec(theta_max_end,runargs,gnostics)\n",
    "\n",
    "nwave = inwavenum.size\n",
    "cfunc = np.reshape(cfunc,[cfunc.shape[1],cfunc.shape[2]])\n",
    "fwhm = 0.005\n",
    "wlen = shiftspec.shape[1]\n",
    "wint =  shiftspec[0,0] - shiftspec[0,wlen-1]\n",
    "# convolve with instrumental profile\n",
    "# start by setting up kernel\n",
    "# First step is finding the array index length of the FWHM\n",
    "disp = wint / wlen\n",
    "gwidth = int((((fwhm / disp) // 2) * 2) +1)\n",
    "# needs to be odd\n",
    "# now get the kernel and convolve\n",
    "gauss = Gaussian1DKernel(gwidth)\n",
    "\n",
    "for ilayer in range (0,press.size):\n",
    "    cfunc[:,ilayer] = convolve(cfunc[:,ilayer],gauss,boundary='extend')\n",
    "\n",
    "tau1_cl_Press = convolve(clphotspec[0],gauss,boundary='extend')[::-1]\n",
    "tau1_oth_Press = convolve(ophotspec[0],gauss,boundary='extend')[::-1]\n",
    "\n",
    "wavenew = shiftspec[0,::-1]\n",
    "press = press.reshape(64,)\n",
    "normfunc = np.zeros_like(cfunc)\n",
    "for iwave in range(0,nwave):\n",
    "    totcont = np.sum(cfunc[iwave,:])\n",
    "    normfunc[iwave,:] = cfunc[iwave,:] / totcont\n",
    "\n",
    "plt.rc('font',family='Times New Roman')\n",
    "fig=plt.figure(dpi=120)\n",
    "plt.rc('font',family='Times New Roman')\n",
    "fig=plt.figure(dpi=320)\n",
    "plt.axis([1.0,15.0,1000,1e-5])\n",
    "\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_xscale('log')\n",
    "#major_ticks = np.arange(1.0,15.,1.0)\n",
    "#minor_ticks = np.arange(1.0,15.,0.5)\n",
    "#ax.set_xticks(major_ticks)                                                       \n",
    "#ax.set_xticks(minor_ticks, minor=True)                                           \n",
    "plt.pcolormesh(wavenew,press,(normfunc[::-1,:].transpose()),cmap='Greys',norm=colors.SymLogNorm(linthresh=0.001,linscale=0.00001,\n",
    "                                              vmin=0., vmax=np.amax(normfunc)))\n",
    "#ax.plot([0.5,20.],[10.**cloud_Top[1],10.**cloud_Top[1]],'c-')\n",
    "plt.plot(wavenew,(tau1_cl_Press),'m-')\n",
    "plt.plot(wavenew,(tau1_oth_Press),'c-')\n",
    "\n",
    "cbar = plt.colorbar(orientation='vertical',norm=colors.Normalize(clip=False),ticks=[1e-3,1e-2,0.1])\n",
    "cbar.ax.set_yticklabels(['<0.1%', '1%', '10%'])\n",
    "cbar.set_label('% of total', rotation=270)\n",
    "plt.ylabel(\"log (P / bar)\")\n",
    "plt.xlabel('Wavelength / Microns')\n",
    "#plt.savefig(runname+'_contribution.png',format='png', dpi=120)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this plots the gas fractions against a solar metallicity chemical equilibrium prediction\n",
    "# IT WILL NEED HACKING FOR YOUR CASE. \n",
    "\n",
    "intab = np.loadtxt(\"2015_06_1060grid_feh_00_co_10.txt\",skiprows=1)\n",
    "test = np.array(intab)\n",
    "test2 = test.reshape(60,18,36)\n",
    "Pgrid = test2[20:21,:,1]\n",
    "Pgrid = Pgrid.reshape(18)\n",
    "Tgrid =  test2[:,10:11,0].reshape(60)\n",
    "abunds= test2[:,:,2:]\n",
    "logP, profT = np.loadtxt(runname+\"_profile.dat\",unpack=True)\n",
    "nlayers = logP.size\n",
    "npress = 18\n",
    "ntemp = 60\n",
    "ngas = 34\n",
    "ab_myP = np.empty([ntemp,nlayers,ngas])\n",
    "for gas in range (0,ngas):\n",
    "    for i in range (0,ntemp):\n",
    "            pfit = interp1d(Pgrid,np.log10(abunds[i,:,gas]))\n",
    "            ab_myP[i,:,gas] = pfit(logP)\n",
    "\n",
    "            \n",
    "# Now get the abundances for this profile....\n",
    "ab = np.empty([nlayers,ngas])\n",
    "for gas in range(0,ngas):\n",
    "    for i in range(0,nlayers):\n",
    "        tfit = InterpolatedUnivariateSpline(Tgrid,ab_myP[:,i,gas],k=1)\n",
    "        ab[i,gas] = tfit(profT[i])\n",
    "vmr = gasVMR\n",
    "plt.rc('font',family='Times New Roman')\n",
    "fig=plt.figure(dpi=320)\n",
    "xl = 1.4\n",
    "yl= 3.8\n",
    "lsize = 10\n",
    "# row and column sharing\n",
    "fig, ((ax1, ax2),(ax3,ax4), (ax5,ax6)) = plt.subplots(3,2, sharex=False, sharey=True)\n",
    "#plt.tight_layout(pad=2.0, w_pad=2.0, h_pad=3.0)\n",
    "plt.subplots_adjust(hspace=0.3)\n",
    "\n",
    "\n",
    "ax1.axis([-12.0,0.,3.0,-4.0])\n",
    "ax1.plot(ab[:,9],logP,'--',color='cyan',linewidth=2)\n",
    "c1, = ax1.plot([vmr[0,1],vmr[0,1]],[2.4,-4.0],'-',color='cyan',label='H2O')\n",
    "ax1.fill_betweenx([2.4,-4.0],[vmr[0,0],vmr[0,0]],[vmr[0,2],vmr[0,2]], facecolor='cyan', alpha=0.1,linewidth=0)\n",
    "\n",
    "ax1.plot(ab[:,11],logP,'--',color='magenta',linewidth=2)\n",
    "c2, = ax1.plot([vmr[1,1],vmr[1,1]],[2.4,-4.0],'-',color='magenta',label='CO')\n",
    "ax1.fill_betweenx([2.4,-4.0],[vmr[1,0],vmr[1,0]],[vmr[1,2],vmr[1,2]], facecolor='magenta', alpha=0.1,linewidth=0)\n",
    "\n",
    "ax1.plot(ab[:,25],logP,'--',color='red',linewidth=2)\n",
    "c3, = ax1.plot([vmr[2,1],vmr[2,1]],[2.4,-4.0],'-',color='red',label='CO2')\n",
    "ax1.fill_betweenx([2.4,-4.0],[vmr[2,0],vmr[2,0]],[vmr[2,2],vmr[2,2]], facecolor='red', alpha=0.1,linewidth=0)\n",
    "\n",
    "ax1.plot(ab[:,10],logP,'--',color='green',linewidth=2)\n",
    "c4, = ax1.plot([vmr[3,1],vmr[3,1]],[2.4,-4.0],'-',color='green',label='CH4')\n",
    "ax1.fill_betweenx([2.4,-4.0],[vmr[3,0],vmr[3,0]],[vmr[3,2],vmr[3,2]], facecolor='green', alpha=0.1,linewidth=0)\n",
    "lgd = ax1.legend(handles=[c1,c2,c3,c4],prop={'size':lsize},loc='best', bbox_to_anchor=(-0.2, 1.0),\n",
    "          ncol=2)\n",
    "\n",
    "      \n",
    "ax2.axis([-12.0,-4.0,3.0,-4.0])      \n",
    "ax2.plot(ab[:,16],logP,'--',color='red',linewidth=2)\n",
    "c3, = ax2.plot([vmr[4,1],vmr[4,1]],[2.4,-4.0],'-',color='red',label='TiO')\n",
    "ax2.fill_betweenx([2.4,-4.0],[vmr[4,0],vmr[4,0]],[vmr[4,2],vmr[4,2]], facecolor='red', alpha=0.1,linewidth=0)\n",
    "ax2.legend(handles=[c3],prop={'size':lsize})\n",
    "\n",
    "ax3.axis([-12.0,-4.0,3.0,-4.0])      \n",
    "ax3.plot(ab[:,17],logP,'--',color='red',linewidth=2)\n",
    "c4, = ax3.plot([vmr[5,1],vmr[5,1]],[2.4,-4.0],'-',color='red',label='VO')\n",
    "ax3.fill_betweenx([2.4,-4.0],[vmr[5,0],vmr[5,0]],[vmr[5,2],vmr[5,2]], facecolor='red', alpha=0.1,linewidth=0)\n",
    "ax3.legend(handles=[c4],prop={'size':lsize})\n",
    "\n",
    "\n",
    "ax4.axis([-12.0,-4.0,3.0,-4.0])      \n",
    "ax4.plot(ab[:,20],logP,'--',color='red',linewidth=2)\n",
    "c5, = ax4.plot([vmr[6,1],vmr[6,1]],[2.4,-4.0],'-',color='red',label='CrH')\n",
    "ax4.fill_betweenx([2.4,-4.0],[vmr[6,0],vmr[6,0]],[vmr[6,2],vmr[6,2]], facecolor='red', alpha=0.1,linewidth=0)\n",
    "ax4.legend(handles=[c5],prop={'size':lsize})\n",
    "\n",
    "ax5.axis([-12.0,-4.0,3.0,-4.0])      \n",
    "ax5.plot(ab[:,19],logP,'--',color='red',linewidth=2)\n",
    "c6, = ax5.plot([vmr[7,1],vmr[7,1]],[2.4,-4.0],'-',color='red',label='FeH')\n",
    "ax5.fill_betweenx([2.4,-4.0],[vmr[7,0],vmr[7,0]],[vmr[7,2],vmr[7,2]], facecolor='red', alpha=0.1,linewidth=0)\n",
    "ax5.legend(handles=[c6],prop={'size':lsize})\n",
    "\n",
    "alks = np.log10((10.**ab[:,22]) + (10.**ab[:,21]))\n",
    "\n",
    "ax6.axis([-10.,-4.0,3.0,-4.0])      \n",
    "ax6.plot(alks,logP,'--',color='red',linewidth=2)\n",
    "c7, = ax6.plot([vmr[8,1],vmr[8,1]],[2.4,-4.0],'-',color='red',label='Na+K')\n",
    "ax6.fill_betweenx([2.4,-4.0],[vmr[8,0],vmr[8,0]],[vmr[8,2],vmr[8,2]], facecolor='red', alpha=0.1,linewidth=0)\n",
    "ax6.legend(handles=[c7],prop={'size':lsize})\n",
    "\n",
    "txt = fig.text(0.5, 0.04, '$\\log f_{gas}$', ha='center')\n",
    "fig.text(0.04, 0.5, '$\\log P / bars$', va='center', rotation='vertical')\n",
    "\n",
    "\n",
    "plt.savefig(runname+'_abundances.png',format='png', dpi=320, bbox_extra_artists=(lgd,txt,),bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}