#!/bin/csh -f

# choose the bands you want (nm)

set b2 = 800
set b1 = 650
set nSamples = 10 # n samples for function
set nReplicates = 10 # how many random samples

#
# set which parameter to fix
#
set relationParam = (1)
# ie LAI here

#
# set any fixed values
#
#set fix = (2 3 4 5 6 7 8 9)
#set fixValue = (0.2 1.0 0.0 1.0 0.005 0.0 40 0.012)

set fix = (9)
set fixValue = (0.012)


# i.e. fix parameter 2 (soil brightness) to be 0.2 here
# i.e. fix parameter 3 (clumping) to be 1.0 here
# i.e. fix parameter 4 (leafEccentricity) to be 0.0 here
# i.e. fix parameter 5 (NLeafLayers) to be 1.0 here
# i.e. fix parameter 6 (DryMatter) to be 0.005 here
# i.e. fix parameter 7 (senescenceProportion) to be 0.0 here
# i.e. fix parameter 8 (chlorophyll) to be 40.0 here
# i.e. fix parameter 9 (leafWater) to be 0.012 here

# some viewing/illumination parameters
set vza = 0
set sza = 45
set saa = 0

# SAVI term
set L = 0.0
# set VI = "RATIO"
set VI = "NDVI"
# 
# names of parameters
#
set name = (LAI soilBrightness clumping leafEccentricity NLeafLayers DryMatter senescenceProportion chlorophyll leafWater)


###################################################################

# vegetation index analysis using kuusk

# set shell variables to find codes
set ARCH = `uname -m`
set path = (. /home/plewis/bin/$ARCH  /home/plewis/bin/csh $path)


# check for required files
set required = (defaults.dat)
foreach i ($required)
  if ( ! -e $i ) then
    echo "required file $i" not found
    exit 0
  endif
end

# use: 
# -N number of leaf layers
# -l wavelength_start wavelength_end wavelength_step

# get data ranges
set lower = (0.001 0.05 0 -1  1 `gawk < defaults.dat '{print $4}'`)
set upper = (10    0.4  2  1  2.5 `gawk < defaults.dat '{print $5}'`)
set range = (`echo $lower $upper | gawk '{for(i=1;i<=NF/2;i++)print $(i+NF/2)-$i;}'`)


mkdir -p data
# define wavelengths for simulations
set w = (`echo $b1 $b2| gawk '{for(i=1;i<=NF;i++)print $i}' | sort -n| gawk '{b[NR]=$1;} END{print b[1],b[2],b[2]-b[1]}'`)

@ rp = 1
set out = data/vi
while ( $rp <= $#relationParam )
set out = "$out.$name[$rp]"
@ rp++
end
set out = "$out.$b1.$b2"
\rm -f $out

set theseparameters = ($relationParam)
@ count = 1
while ( $count <= $nSamples )
  @ rp = 1
  echo -n "parameter: "
  while ( $rp <= $#relationParam )
    set which = $relationParam[$rp]
    set theseparameters[$rp] = `echo $lower[$which] $upper[$which] | gawk '{d=($2-$1)/(1.*n);print $1 + d*(i-1.0);}' n=$nSamples i=$count`
    echo -n "P${which}: $name[$which] = $theseparameters[$rp] "
    @ rp++
  end
  echo ": bands: $w"
  @ n = 1
  set params = ($lower)
  while ( $n <= $nReplicates)
    echo -n "...$n"
    # randomise variables
    set rand = (`gawk -v np=$#name  'BEGIN{srand();for(i=1;i<=np;i++)print rand()}'`)
    set params = ($rand)
    @ i = 1
    while ( $i <= $#params )
      set params[$i] = `echo $lower[$i] $upper[$i] $rand[$i] | gawk '{print $1+$3*($2-$1)}'`
      @ i++
    end
    # set the control variables
    @ rp = 1
    set P = 1.0
    while ( $rp <= $#relationParam )
      set params[$relationParam] = $theseparameters[$rp]
      set P = `echo "$P*$theseparameters[$rp]" | \bc -l`
      @ rp++
    end
    
    if ( $#fix ) then
        @ k = 1
        while ( $k <= $#fix )
            set this = $fix[$k]
            set that = $fixValue[$k]
            set params[$this] = $that
            @ k++
        end
    endif

    set test = `echo $params[4] |gawk '{if($1<0)print 0;else print 1;}'`
    set testp = `echo $params[4] |gawk '{if($1<0)print -$1;else print $1;}'`
    if ( $test == 0 ) then
	set ee = "-thm 90 -ee $testp"
    else
	set ee = "-thm 0 -ee $testp"
    endif
    BRDF.ip -view $vza $vza 1 -solar_azimuth $sza -solar_zenith $sza -wavelength $w | kuusk -L $params[1] $ee -cAb $params[8] -rsl1 $params[2] -cW $params[9] -N $params[5] -clz $params[3] -cP $params[6] -cC $params[7] -data /home/plewis/src/kuusk/pigments.dat > tmp.$$

    set omega = `head -2 < tmp.$$ | tail -1 | gawk '{for(i=5;i<=NF;i++)print $i}'`
    \rm -f tmp.$$
    # pull out refl for 2 wavelengths of interest
    if ( $VI == "NDVI" ) then
      set vi = `echo $omega | gawk '{print (1.0+L)*($2-$1)/($2+$1+L)}' L=$L`
    else
      # ratio
       set vi = `echo $omega | gawk '{print ($2-L)/$1}' L=$L`
    endif     
    echo $P $vi $omega >> $out

    @ n++
  end
  echo ""
  @ count++
end
sort -n < $out > tmp.$$
mv tmp.$$ $out.$VI.L$L
echo "plot of VI against parameter $relationParam ($name[$relationParam])"
generate_graph $out.$VI.L$L -noline -xlabel $name[$relationParam] -ylabel VI

echo "scatter plot of reflectance at $b2 nm against $b1 nm"

generate_graph $out.$VI.L$L -noline -xlabel "refletance at $b1 nm"  -ylabel "refletance at $b2 nm" -u 3:4