~~~python
import numpy as np
import matplotlib.pyplot as plt
import htrdrPy as htrdr
import scipy.constants as cst
import json
~~~

## Read data and create Data objects

```python
radius = 2.575e6 # m

data = np.load("data/1D_spetral_data.npz")
data = dict(data)
for key, array in data.items():
    print(key, ", shape : ", array.shape)

nWavelength, nLevel, nCoeff = data["absorption (m-1)"].shape
nAngle = len(data["angles (°)"])
wavelengths = data["wavelength"]

data.update({
        'nLevel' : nLevel,
        'nWavelength' : nWavelength,
        'nCoeff' : nCoeff,
        'nAngle' : nAngle,
        })

print(data.keys())

dataSp= htrdr.Data(radius = radius, nTheta=80, nPhi=80, name="Sphere")

# In plan parallel, the radius is used to define the x and y extension of the
# atmosphere. Make sure this value is large enough
dataPP = htrdr.Data(radius = 1e9, name="PP")

dataSp.makeMixture(data, dim=1)
dataPP.makeMixture(data, dim=0)

brdf = {
        "kind": 'lambertian',
        "albedo": data['surface albedo'],
        "wavelengths": wavelengths
        }
dataPP.makeGroundFrom1D_PP(91, brdf)
dataSp.makeGroundFrom1D(91, brdf)
```


## Write input files and VTK files

```python
dataPP.writeInputs()
dataSp.writeInputs()

# VTK files are not necessary and only serves to visualize the data that will be
# send to htrdr-planets
dataPP.writeVTKfiles()
dataSp.writeVTKfiles()
```

## Calculate the observation geometry

```python
with open("data/observation.json", 'r') as f: obs = json.loads(f.read())
cameraDist = 1000 * cst.kilo    # m
solDist = 10 * cst.au   # m
solRad = 7e8    # m
solTemp = 5800  # K

geomPP = htrdr.Geometry(case='PP')
geomPP.setImage([64,64], 10)
geomPP.makeGeomFromAPIE(obs, 0, cameraDist, solDist, solRad, srcTemp=solTemp)
geomPP.exportGeometry()

geomSp = htrdr.Geometry(case='Sphere')
geomSp.setImage([64,64], 10)
geomSp.makeGeomFromAPIE(obs, radius, cameraDist, solDist, solRad, srcTemp=solTemp)
geomSp.exportGeometry()
```

## Generate the scripts

```python
scriptPP = htrdr.Script(case='Spectrum',
        MPIcmd="mpirun -np 4 --map-by socket:PE=9")
scriptPP.reflectanceSpectrum(geomPP, "sw", wavelengths)

scriptSp = htrdr.Script(case='Spectrum',
        MPIcmd="mpirun -np 4 --map-by socket:PE=9")
scriptSp.reflectanceSpectrum(geomSp, "sw", wavelengths)
```

## Start the scripts

```python
scriptPP(dataPP)
scriptSp(dataSp)
```

## Start the post-process

```python
htrdr.Postprocess(scriptPP)
htrdr.Postprocess(scriptSp)
```

## Plot results

```python
with open("results_PP/reflectance_spectrum_Spectrum.json", 'r') as f:
    resPP = json.loads(f.read())

with open("results_Sphere/reflectance_spectrum_Spectrum.json", 'r') as f:
    resSp= json.loads(f.read())

with open("data/spectrum.txt", 'r') as f:
    f.readline()
    wvl = []
    IF_obs   = []
    IF_SHDOM = []
    IF_MCC   = []
    for line in f:
        l = line.split()
        wvl     .append(l[0])
        IF_obs  .append(l[2])
        IF_SHDOM.append(l[3])
        IF_MCC  .append(l[4])

wvl      = np.array(wvl     , dtype=float)
IF_obs   = np.array(IF_obs  , dtype=float)
IF_SHDOM = np.array(IF_SHDOM, dtype=float)
IF_MCC   = np.array(IF_MCC  , dtype=float)

print(wvl)
print(IF_obs)

fig, ax = plt.subplots()

ax.plot(wvl, IF_obs  , label="Observed", ls="", marker="o")
ax.plot(wvl, IF_SHDOM, label="SHDOMPP", ls="", marker="o")
ax.plot(wvl, IF_MCC  , label="MCCSPHE", ls="", marker="o")

ax.errorbar(np.array(resPP['wavelength'])/cst.micro,
            np.array(resPP['reflectance spectrum']),
            yerr = 3 * np.array(resPP['reflectance std deviation']),
            label = "plan-parallel")

ax.errorbar(np.array(resSp['wavelength'])/cst.micro,
            np.array(resSp['reflectance spectrum']),
            yerr = 3 * np.array(resSp['reflectance std deviation']),
            label = "spherical")

ax.set_yscale('log')
ax.set_ylabel('Reflectance')
ax.set_xlabel('Wavlength [µm]')

ax.legend()

fig.savefig('Spectrum.png')
```
