diff --git a/README.md b/README.md
new file mode 100644
index 0000000..480f067
--- /dev/null
+++ b/README.md
@@ -0,0 +1,103 @@
+
+
+
+
+
+
+
+---
+
+# PROSAIL Python Bindings
+
+#### J Gomez-Dans (NCEO & UCL) ``j.gomez-dans@ucl.ac.uk``
+
+
+## Description
+
+This repository contains the Python bindings to the PROSPECT and SAIL leaf and
+canopy reflectance models. The code is written in FORTRAN. The original fortran
+code was downloaded from [Jussieu](http://teledetection.ipgp.jussieu.fr/prosail/).
+I only added a function to simplify writing the wrappers, and you should go to
+that page to get newer versions of the code. A recent review of the use of both
+models is availabe in [this paper](http://webdocs.dow.wur.nl/internet/grs/Workshops/Environmental_Applications_Imaging_Spectroscopy/12_Jacquemoud_Prospect/IEEE_Jacquemoud_PROSPECT.pdf)_.
+
+
+## Installing the bindings
+
+The installation of the bindings is quite straightforward: unpack the distribution
+and run the following command
+
+ python setup.py install
+
+You can usually install it to your user's directory (if you haven't got superuser
+privileges) by
+
+ python setup.py install --user
+
+#### **Note**
+
+
+You will need a working FORTRAN compiler. I have only tested this with GCC on Linux, but it should work on other systems. You can also pass optimisation flags to the compiler:
+
+ python setup.py config_fc --fcompiler=gnu95 --arch=-march=native --opt=-O3 install --user
+
+## Using the bindings
+
+The bindings offer several functions, which will be described in detail below:.
+
+* ``run_prospect``: This function runs the PROSPECT 5B model in Feret et al 2008. The input parameters are the usual ``(n,cab,car,cbrown,cw,cm)`` (e.g. leaf layers, leaf Chlorophyll concentration, leaf Carotenoid concentration, leaf senescent fraction, Equivalent leaf water, leaf dry matter). It returns a spectrum between 400 and 2500 nm.
+* ``run_sail``: The SAILh model, which in this case requires leaf reflectance and transmittance to be fed to the model (e.g. you have already measured these spectra in the field). The rest of the parameters are ``(refl,trans,lai,lidfa,lidfb,rsoil,psoil,hspot,tts,tto,psi,typelidf)``. Additionally, there are two optional parameters, ``soil_spectrum1``, ``soil_spectrum2``, which allow you to set the soil spectra (otherwise, some default spectra get used). Output is a spectrum in the 400-2500 nm range.
+* ``run_prosail``: PROSPECT_5B and SAILh coupled together, with input parameters given by ``(n,cab,car,cbrown,cw,cm,lai,lidfa,lidfb,rsoil,psoil,hspot,tts,tto,psi,typelidf)``. Output is a spectrum in the 400-2500 nm range.
+
+
+## The parameters
+
+The parameters used by the models and their units are introduced below:
+
+| Parameter | Description of parameter | Units |Typical min | Typical max |
+|-------------|---------------------------------|--------------|------------|-------------|
+| N | Leaf structure parameter | N/A | 0.8 | 2.5 |
+| cab | Chlorophyll a+b concentration | ug/cm2 | 0 | 80 |
+| caw | Equivalent water thickiness | cm | 0 | 200 |
+| car | Carotenoid concentration | ug/cm2 | 0 | 20 |
+| cbrown | Brown pigment | NA | 0 | 1 |
+| cm | Dry matter content | g/cm2 | 0 | 200 |
+| lai | Leaf Area Index | N/A | 0 | 10 |
+| lidfa | Leaf angle distribution | N/A | - | - |
+| lidfb | Leaf angle distribution | N/A | - | - |
+| psoil | Dry/Wet soil factor | N/A | 0 | 1 |
+| rsoil | Soil brigthness factor | N/A | - | - |
+| hspot | Hotspot parameter | N/A | - | - |
+| tts | Solar zenith angle | deg | 0 | 90 |
+| tto | Observer zenith angle | deg | 0 | 90 |
+| phi | Relative azimuth angle | deg | 0 | 360 |
+| typelidf | Leaf angle distribution type | Integer | - | - |
+
+### Specifying the leaf angle distribution
+
+The parameter ``typelidf`` regulates the leaf angle distribution family being used. The following options are understood:
+
+* ``typelidf = 1``: use the two parameter LAD parameterisation, where ``a`` and ``b`` control the average leaf slope and the distribution bimodality, respectively. Typical distributions
+are given by the following parameter choices:
+
+| LIDF type | ``LIDFa`` | ``LIDFb`` |
+|--------------|-----------|------------------|
+| Planophile | 1 | 0 |
+| Erectophile| -1 | 0 |
+| Plagiophile| 0 | -1 |
+| Extremophile| 0 | 1 |
+| Spherical | -0.35 | -0.15 |
+| Uniform | 0 | 0 |
+
+* ``typelidf = 2`` Ellipsoidal distribution, where ``LIDFa`` parameter stands for mean leaf angle (0 degrees is planophile, 90 degrees is erectophile). ``LIDFb`` parameter is ignored.
+
+### The soil model
+
+The soil model is a fairly simple linear mixture model, where two spectra are mixed and then a brightness term added:
+
+ rho_soil = rsoil*(psoil*soil_spectrum1+(1-psoil)*soil_spectrum2)
+
+
+The idea is that one of the spectra is a dry soil and the other a wet soil, so soil moisture is then contorlled by ``psoil``. ``rsoil`` is just a brightness scaling term.
+
+
diff --git a/README.rst b/README.rst
deleted file mode 100644
index 4fc8c7f..0000000
--- a/README.rst
+++ /dev/null
@@ -1,119 +0,0 @@
-PROSAIL
-==========
-
-Description
---------------
-
-This repository contains the Python bindings to the PROSPECT and SAIL leaf and
-canopy reflectance models. The code is written in FORTRAN. The original fortran
-code was downloaded from `Jussieu `_.
-I only added a function to simplify writing the wrappers, and you should go to
-that page to get newer versions of the code. A recent review of the use of both
-models is availabe in `this paper `_.
-
-
-Installing the bindings
--------------------------
-
-The installation of the bindings is quite straightforward: unpack the distribution
-and run the following command
-
- python setup.py install
-
-You can usually install it to your user's directory (if you haven't got superuser
-privileges) by
-
- python setup.py install --user
-
-Note
-*******
-
-You will need a working FORTRAN compiler. I have only tested this with GCC on Linux, but it should work on other systems. You can also pass optimisation flags to the compiler:
-
- python setup.py config_fc --fcompiler=gnu95 --arch=-march=native --opt=-O3 install --user
-
-
-Using the bindings
----------------------
-
-The bindings offer a single function, ``run_prosail``.
-
- retval = run_prosail(n,cab,car,cbrown,cw,cm,lai,lidfa,lidfb,psoil,hspot,tts,tto,psi)
-
-The parameters that are fed into ``run_prosail`` are
-
-+-------------+---------------------------------+--------------+------------+-------------+
-| Parameter | Description of parameter | Units |Typical min | Typical max |
-+=============+=================================+==============+============+=============+
-| N | Leaf structure parameter | N/A | 0.8 | 2.5 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| cab | Chlorophyll a+b concentration | ug/cm2 | 0 | 200 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| caw | Equivalent water thickiness | cm | 0 | 200 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| car | Carotenoid concentration | ug/cm2 | 0 | 200 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| cbrown | Brown pigment | NA | 0 | 1 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| cm | Dry matter content | g/cm2 | 0 | 200 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| lai | Leaf Area Index | N/A | 0.1 | 10 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| lidfa | Leaf angle distribution | N/A | - | - |
-+-------------+---------------------------------+--------------+------------+-------------+
-| lidfb | Leaf angle distribution | N/A | - | - |
-+-------------+---------------------------------+--------------+------------+-------------+
-| psoil | Dry/Wet soil factor | N/A | 0 | 1 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| hspot | Hotspot parameter | N/A | 0 | 0. |
-+-------------+---------------------------------+--------------+------------+-------------+
-| tts | Solar zenith angle | deg | 0 | 90 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| tto | Observer zenith angle | deg | 0 | 90 |
-+-------------+---------------------------------+--------------+------------+-------------+
-| phi | Relative azimuth angle | deg | 0 | 360 |
-+-------------+---------------------------------+--------------+------------+-------------+
-
-``lidfa`` and ``lidfb`` parameters control the leaf angle distribution. Typical distributions
-are given by the following parameter choices:
-
-+--------------+-----------+------------------+
-|LIDF type | lidfa | lidfb |
-+==============+===========+==================+
-|Planophile | 1 | b |
-+--------------+-----------+------------------+
-| Erectophile| -1 | 0 |
-+--------------+-----------+------------------+
-| Plagiophile| 0 | -1 |
-+--------------+-----------+------------------+
-| Extremophile| 0 | 1 |
-+--------------+-----------+------------------+
-| Spherical | -0.35 | -0.15 |
-+--------------+-----------+------------------+
-| Uniform | 0 | 0 |
-+--------------+-----------+------------------+
-
-
-
-
-A simple example of using bindings would be
-
- cm=0.009
- cab = 80
- car = 15
- cbrown = 0
- cw = 0.001
- lidfa = 0
- lidfb = 0
- psoil = 0
- hspot = 0.01
- tts = 30
- tto = 10
- phi = 0
- lai = np.arange ( 0, 5, 0.2 )
- for l in lai:
- plt.plot ( run_prosail(1.5,cab,car,cbrown,cw,cm,l,lidfa,lidfb,psoil,hspot,tts,tto,psi))
-
-This results in a simulation of surface reflectance (from 400 to 2500 nm) as a function of LAI and under the other parameters' prescribed values. This yields the following
-
-.. image:: http://i.imgur.com/2Hh0z.png
diff --git a/prosail/__init__.py b/prosail/__init__.py
index 346dd8b..51f48c9 100644
--- a/prosail/__init__.py
+++ b/prosail/__init__.py
@@ -1,5 +1,141 @@
import numpy as np
-from prosail_fortran import run_sail, run_prosail, prospect_5b
+from prosail_fortran import run_sail as sail, run_prosail as prosail
+from prosail_fortran import prospect_5b
+from prosail_fortran import mod_dataspec_p5b as spectral_libs
+
+def run_prosail (n,cab,car,cbrown,cw,cm,lai,lidfa,lidfb,rsoil,psoil,hspot,
+ tts,tto,psi,typelidf,
+ soil_spectrum1=None,soil_spectrum2=None ):
+ """Run the PROSPECT_5B and SAILh radiative transfer models. The soil
+ model is a linear mixture model, where two spectra are combined together as
+
+ rho_soil = rsoil*(psoil*soil_spectrum1+(1-psoil)*soil_spectrum2)
+ By default, ``soil_spectrum1`` is a dry soil, and ``soil_spectrum2`` is a
+ wet soil, so in that case, ``psoil`` is a surface soil moisture parameter.
+ ``rsoil`` is a soil brightness term. You can provide one or the two
+ soil spectra if you want. The soil spectra must be defined
+ between 400 and 2500 nm with 1nm spacing.
+
+ Parameters
+ ----------
+ n: float
+ Leaf layers
+ cab: float
+ leaf chlorophyll concentration
+ car: float
+ leaf carotenoid concentration
+ cbrown: float
+ senescent pigment
+ cw: float
+ equivalent leaf water
+ cm: float
+ leaf dry matter
+ lai: float
+ leaf area index
+ lidfa: float
+ a parameter for leaf angle distribution. If ``typliedf``=2, average
+ leaf inclination angle.
+ lidfb: float
+ b parameter for leaf angle distribution. If ``typelidf``=2, ignored
+ rsoil: float
+ Soil scalar
+ psoil: float
+ Soil scalar
+ tts: float
+ Solar zenith angle
+ tto: float
+ Sensor zenith angle
+ psi: float
+ Relative sensor-solar azimuth angle ( saa - vaa )
+ soil_spectrum1: 2101-element array
+ First component of the soil spectrum
+ soil_spectrum2: 2101-element array
+ Second component of the soil spectrum
+
+ Returns
+ --------
+ Directional surface reflectance between 400 and 2500 nm
+
+
+ """
+
+ if soil_spectrum1 is not None:
+ assert ( len(soil_spectrum1) == 2101 )
+ else:
+ soil_spectrum1 = spectral_libs.rsoil1
+
+ if soil_spectrum2 is not None:
+ assert ( len(soil_spectrum1) == 2101 )
+ else:
+ soil_spectrum2 = spectral_libs.rsoil1
+
+ rho = prosail (n,cab,car,cbrown,cw,cm,lai,lidfa,lidfb,rsoil,psoil,hspot,
+ tts,tto,psi,typelidf, soil_spectrum1, soil_spectrum2 )
+ return rho
+
+def run_sail (refl,trans,lai,lidfa,lidfb,rsoil,psoil,hspot,tts,tto,psi,typelidf,
+ soil_spectrum1=None,soil_spectrum2=None ):
+ """Run the SAILh radiative transfer model. The soil model is a linear
+ mixture model, where two spectra are combined together as
+
+ rho_soil = rsoil*(psoil*soil_spectrum1+(1-psoil)*soil_spectrum2)
+
+ By default, ``soil_spectrum1`` is a dry soil, and ``soil_spectrum2`` is a
+ wet soil, so in that case, ``psoil`` is a surface soil moisture parameter.
+ ``rsoil`` is a soil brightness term. You can provide one or the two
+ soil spectra if you want. The soil spectra, and leaf spectra must be defined
+ between 400 and 2500 nm with 1nm spacing.
+
+ Parameters
+ ----------
+ refl: 2101-element array
+ Leaf reflectance
+ trans: 2101-element array
+ leaf transmittance
+ lai: float
+ leaf area index
+ lidfa: float
+ a parameter for leaf angle distribution. If ``typliedf``=2, average
+ leaf inclination angle.
+ lidfb: float
+ b parameter for leaf angle distribution. If ``typelidf``=2, ignored
+ rsoil: float
+ Soil scalar
+ psoil: float
+ Soil scalar
+ tts: float
+ Solar zenith angle
+ tto: float
+ Sensor zenith angle
+ psi: float
+ Relative sensor-solar azimuth angle ( saa - vaa )
+ soil_spectrum1: 2101-element array
+ First component of the soil spectrum
+ soil_spectrum2: 2101-element array
+ Second component of the soil spectrum
+
+ Returns
+ --------
+ Directional surface reflectance between 400 and 2500 nm
+
+
+ """
+
+
+ if soil_spectrum1 is not None:
+ assert ( len(soil_spectrum1) == 2101 )
+ else:
+ soil_spectrum1 = spectral_libs.rsoil1
+
+ if soil_spectrum2 is not None:
+ assert ( len(soil_spectrum1) == 2101 )
+ else:
+ soil_spectrum2 = spectral_libs.rsoil1
+
+ rho = sail (refl,trans,lai,lidfa,lidfb,rsoil,psoil,hspot,tts,tto,psi,typelidf,
+ soil_spectrum1, soil_spectrum2 )
+ return rho
+
def trans_prosail ( N, cab, car, cbrown, cw, cm, lai, lidfa, lidfb, rsoil, psoil, \
diff --git a/prosail/dataSpec_P5B.f90 b/prosail/dataSpec_P5B.f90
index fb2135c..856f949 100644
--- a/prosail/dataSpec_P5B.f90
+++ b/prosail/dataSpec_P5B.f90
@@ -2268,4 +2268,4 @@ MODULE MOD_dataSpec_P5B
4.960E-02,4.960E-02,4.960E-02,4.960E-02,4.959E-02,4.959E-02,4.944E-02,4.930E-02,4.915E-02,4.900E-02,&
4.885E-02/
- END
\ No newline at end of file
+ END
diff --git a/prosail/prosail_fortran.pyf b/prosail/prosail_fortran.pyf
index e336028..b8a5b89 100644
--- a/prosail/prosail_fortran.pyf
+++ b/prosail/prosail_fortran.pyf
@@ -18,7 +18,7 @@ python module prosail_fortran ! in
real*8 dimension(2101) :: es
integer, parameter,optional :: nw=2101
end module mod_dataspec_p5b
- subroutine run_prosail(n,cab,car,cbrown,cw,cm,lai,lidfa,lidfb,rsoil,psoil,hspot,tts,tto,psi,typelidf,retval, soil_spectrum1, soil_spectrum2) ! in :prosail_fortran:run_prosail.f90
+ subroutine run_prosail(n,cab,car,cbrown,cw,cm,lai,lidfa,lidfb,rsoil,psoil,hspot,tts,tto,psi,typelidf,soil_spectrum1, soil_spectrum2, retval) ! in :prosail_fortran:run_prosail.f90
use mod_sail
use mod_angle
use mod_flag_util
@@ -41,11 +41,11 @@ python module prosail_fortran ! in
real*8 intent(in) :: tto
real*8 intent(in) :: psi
integer intent(in) :: typelidf
- real*8, dimension(2101), intent(in), optional :: soil_spectrum1
- real*8, dimension(2101), intent(in), optional :: soil_spectrum2
+ real*8, dimension(2101), intent(in) :: soil_spectrum1
+ real*8, dimension(2101), intent(in) :: soil_spectrum2
real*8 dimension(2101),intent(out) :: retval
end subroutine run_prosail
- subroutine run_sail(refl,trans,lai,lidfa,lidfb,rsoil,psoil,hspot,tts,tto,psi,typelidf,retval) ! in :prosail_fortran:run_prosail.f90
+ subroutine run_sail(refl,trans,lai,lidfa,lidfb,rsoil,psoil,hspot,tts,tto,psi,typelidf,soil_spectrum1, soil_spectrum2, retval) ! in :prosail_fortran:run_prosail.f90
use mod_sail
use mod_angle
use mod_flag_util
@@ -64,6 +64,8 @@ python module prosail_fortran ! in
real*8 intent(in) :: tto
real*8 intent(in) :: psi
integer intent(in) :: typelidf
+ real*8,dimension(2101),intent(in) :: soil_spectrum1
+ real*8,dimension(2101),intent(in) :: soil_spectrum2
real*8 dimension(2101),intent(out) :: retval
end subroutine run_sail
subroutine prospect_5b(n,cab,car,cbrown,cw,cm,rt) ! in :prosail_fortran:prospect_5B.f90
diff --git a/prosail/run_prosail.f90 b/prosail/run_prosail.f90
index 5aa8aef..9cba7c9 100644
--- a/prosail/run_prosail.f90
+++ b/prosail/run_prosail.f90
@@ -1,7 +1,7 @@
SUBROUTINE run_prosail ( N, Cab, Car, Cbrown, Cw, Cm, lai, LIDFa, LIDFb, &
- rsoil, psoil, hspot, tts, tto, psi, TypeLidf, retval, &
- soil_spectrum1, soil_spectrum2 )
+ rsoil, psoil, hspot, tts, tto, psi, TypeLidf, &
+ soil_spectrum1, soil_spectrum2, retval )
USE MOD_ANGLE ! defines pi & rad conversion
USE MOD_staticvar ! static variables kept in memory for optimization
@@ -21,8 +21,9 @@ SUBROUTINE run_prosail ( N, Cab, Car, Cbrown, Cw, Cm, lai, LIDFa, LIDFb, &
INTEGER, intent(in) :: TypeLidf
REAL*8,ALLOCATABLE,SAVE :: resh(:),resv(:)
REAL*8,ALLOCATABLE,SAVE :: rsoil0(:),PARdiro(:),PARdifo(:)
- REAL*8, dimension(nw), intent(in), optional :: soil_spectrum1
- REAL*8, dimension(nw), intent(in), optional :: soil_spectrum2
+ REAL*8, dimension(nw), intent(in) :: soil_spectrum1
+ REAL*8, dimension(nw), intent(in) :: soil_spectrum2
+
INTEGER :: ii
REAL*8 :: ihot, skyl
! ANGLE CONVERSION
@@ -85,13 +86,11 @@ SUBROUTINE run_prosail ( N, Cab, Car, Cbrown, Cw, Cm, lai, LIDFa, LIDFb, &
ALLOCATE (rsoil0(nw))
!psoil = 1. ! soil factor (psoil=0: wet soil / psoil=1: dry soil)
! rsoil : soil brightness term
- if ( present(soil_spectrum1) ) then
- Rsoil1 = soil_spectrum1
- endif
- if ( present ( soil_spectrum2 ) ) then
- Rsoil2 = soil_spectrum2
- endif
- rsoil0=rsoil*(psoil*Rsoil1+(1-psoil)*Rsoil2)
+
+
+
+ rsoil0=rsoil*(psoil*soil_spectrum1+(1-psoil)*soil_spectrum2)
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 4SAIL canopy structure parm !!
@@ -147,7 +146,8 @@ END subroutine run_prosail
SUBROUTINE run_sail ( refl, trans, lai, LIDFa, LIDFb, &
- rsoil, psoil, hspot, tts, tto, psi, TypeLidf, retval )
+ rsoil, psoil, hspot, tts, tto, psi, TypeLidf, retval, &
+ soil_spectrum1, soil_spectrum2 )
USE MOD_ANGLE ! defines pi & rad conversion
USE MOD_staticvar ! static variables kept in memory for optimization
@@ -167,6 +167,9 @@ SUBROUTINE run_sail ( refl, trans, lai, LIDFa, LIDFb, &
INTEGER, intent(in) :: TypeLidf
REAL*8,ALLOCATABLE,SAVE :: resh(:),resv(:)
REAL*8,ALLOCATABLE,SAVE :: rsoil0(:),PARdiro(:),PARdifo(:)
+ REAL*8, dimension(nw), intent(in) :: soil_spectrum1
+ REAL*8, dimension(nw), intent(in) :: soil_spectrum2
+
INTEGER :: ii
REAL*8 :: ihot, skyl
! ANGLE CONVERSION
@@ -229,7 +232,8 @@ SUBROUTINE run_sail ( refl, trans, lai, LIDFa, LIDFb, &
ALLOCATE (rsoil0(nw))
!psoil = 1. ! soil factor (psoil=0: wet soil / psoil=1: dry soil)
! rsoil : soil brightness term
- rsoil0=rsoil*(psoil*Rsoil1+(1-psoil)*Rsoil2)
+
+ rsoil0=rsoil*(psoil*soil_spectrum1+(1-psoil)*soil_spectrum2)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 4SAIL canopy structure parm !!
@@ -245,6 +249,7 @@ SUBROUTINE run_sail ( refl, trans, lai, LIDFa, LIDFb, &
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! CALL PRO4SAIL !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ print*, refl(400),tau(400),rsoil0(400)
CALL SAIL(refl, trans,LIDFa,LIDFb,TypeLIDF,LAI,hspot,tts,tto,psi,rsoil0)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -282,4 +287,4 @@ SUBROUTINE run_sail ( refl, trans, lai, LIDFa, LIDFb, &
!WRITE( *,'(i4,f10.6)') (lambda(ii),run_prosail(ii), ii=1,nw)
END subroutine run_sail
-
\ No newline at end of file
+