<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 PyDISORT
Step3:  Define methods for generating an atmospheric model of Titan and performing RT calculation
Note the parameters that are specified in the methods below. For example, the k-coefficients file,
the gas compositions, and the surface reflectivity.
Step4:  Parallel (multi-core) execution
Requires ipyparallel and starting clusters in IPython Clusters tab of notebook to make use of multilpe cores.
Step5:  Compare calculations with observed VIMS spectrum
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import os
from time import time
from ipyparallel import Client
os.environ['RTDATAPATH'] = '/Users/mate/g/rt/data/refdata/'
import atmosphere as atm
from atmosphere.rt import pydisort
atm.refdata.setup_directory()
Explanation: PyDISORT: Example calculation of VIMS spectrum
&#109;&#97;&#116;&#101;&#64;&#98;&#101;&#114;&#107;&#101;&#108;&#101;&#121;&#46;&#101;&#100;&#117;
The following notebook gives an example of how to setup the atmospheric 
structure dictionary and use it to calculate a spectrum at the VIMS 
dispersion (wavelength) scale. The pydisort module has been updated to include
a parallel (multiprocessing) implementation of the call to the monochromatic CDISORT
code. The requires ipyparallel. 
The notebook below demonstrates how to change a few of the default parameters,
run the radiative transfer calculation in either single process or multiprocessing
mode and then illustrates the comparison of the model to some VIMS observations.
Specify data directory and download/confirm reference data
End of explanation
def create_atmosphere_model(**kw):
    
    Setup layered atmospheric model using a default physical/thermal structure
    determined by HASI, and composition from the GCMS, with aerosol 
    scattering propertier from DISR.
    
    The titan dictionary contains the values used to determine the opacity
    and scattering albedo of each atmospheric layer as a function of
    wavelength. The dispersion axis is constrained by the k-coefficients that
    are used, as specified by the methane kc_file.
    
    
    ts = time()
    titan = atm.structure.set_HASI_structure(nlev=21,method='split_at_tropopause')
    atm.composition.set_abundances(titan, trace_gas={'m_H2':0.001})
    atm.gas_opacity.set_methane(titan, 
#        kc_file='/Users/mate/g/rt/data/refdata/gas_opacity/kc_CH4.VIMS.v08.fits',
        kc_file='/Users/mate/data/k/coeff/CH4/kc_CH4.HST.v01.fits',                                
        )
    atm.gas_opacity.set_cia(titan)
    atm.aerosol.set_opacity(titan)
    DISR = atm.aerosol.fit_DISR_phases()
    atm.aerosol.set_aerosol_phase_moments(titan, DISR, nmom=32)   
    t_setup = time()-ts
    
    titan['haze'].update({'ssalb':0.96})
    titan.update({'radius':2575.,
                  'rsurf':0.10,
                  't_setup':t_setup,
                 })
    if 'rsurf' in kw: titan.update({'rsurf':kw['rsurf']})
    if 'verbose' in kw and kw['verbose']:
        pstr = 'PyDISORT Titan atmosphere structure and opacity setup: {:6.1f} sec'
        print(pstr.format(titan['t_setup']))    
    
    return titan
def setup_VIMS_calc(**kw):
    Specify the viewing geometry and wavelegth range for the
    radiative transfer calculation.
    titan = create_atmosphere_model(**kw)
    titan.update({'rt':{'spher':False,
                        'wav_range':(2.0,2.40),
                        'view':{'umu0':0.99,'umue': 0.90,'phi0': 10.0,'phie': 11.0},
                       }})
    for k in ['view','wav_range']: 
        if k in kw: titan['rt'].update({k:kw[k]})
            
    fi = lambda array, v: abs(array-v).argmin()
    wav_indices = lambda array, wavs: tuple([fi(array, mu) for mu in wavs])
    wav_mn, wav_mx = wav_indices(titan['wavelength'], titan['rt']['wav_range'])
    titan['rt'].update({'wav_mn':wav_mn,
                        'wav_mx':wav_mx,
                        'nlam':wav_mx-wav_mn+1,
                        })          
    return titan
titan = setup_VIMS_calc(wav_range=(0.8900,0.8901), verbose=True ) 
pydisort.calculate_spectrum(titan)
titan['rt']
def spectest(model):
    for i in range(model['rt']['wav_mn'],model['rt']['wav_mx']):
        print(i, model['wavelength'][i], model['spectrum'][i])
        
spectest(titan)
from astropy.io import fits
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(titan['wavelength'], titan['spectrum'])
ax.set_ylim(0.18,0.20)
Explanation: Define methods for generating an atmospheric model of Titan and performing RT calculation
Note the parameters that are specified in the methods below. For example, the k-coefficients file,
the gas compositions, and the surface reflectivity.
End of explanation
from ipyparallel import parallel, Client
%%px 
import os
os.environ['RTDATAPATH'] = '/Users/mate/g/rt/data/refdata/'
pwd
titan = setup_VIMS_calc(wav_range=(0.8850,0.8890))
pydisort.ipcluster_spectrum_calculation(titan)
titan['rt']
Explanation: Parallel (multi-core) execution
Requires ipyparallel and starting clusters in IPython Clusters tab of notebook to make use of multilpe cores.
End of explanation
VIMS_test = setup_VIMS_calc(
                rsurf=0.25,
                wav_range=(1.5,3.5),
                view={'umu0':np.cos(73.13*(np.pi/180)),
                      'umue':np.cos(51.93*(np.pi/180)),
                      'phi0': 272.4,
                      'phie': 360-111.4},
                verbose=True,
                ) 
pydisort.cluster_spectrum_calculation(VIMS_test)
VIMS_test['rt']
%matplotlib inline
import matplotlib.pyplot as plt
url = 'http://w.astro.berkeley.edu/~madamkov/refdata/test/VIMS.dat'
obs = np.genfromtxt(url,names=['wav','ref','cloud'])
fig, ax = plt.subplots()
ax.plot(obs['wav'], obs['cloud'], 'b-', drawstyle='steps-mid', label='cloud')
ax.plot(obs['wav'], obs['ref'], 'g-', drawstyle='steps-mid', label='reference')
ax.plot(VIMS_test['wavelength'], VIMS_test['spectrum'], 'r-', drawstyle='steps-mid', label='PyDISORT') 
ax.set_xlim(1.2,3.5) ;
ax.legend(loc='best')
ax.set_xlabel('wavelength (um)') ; ax.set_ylabel('I/F') ;
Explanation: Compare calculations with observed VIMS spectrum
End of explanation
<END_TASK>