|
| 1 | +## BASH shell code to calculate the theoretical normal reflectivity of a 3-layer silver scale ultrastructure. The scale layer thickness parameters can be varied one at a time, any two or all three simultaneously. |
| 2 | +## Author: Vinodkumar Saranathan (https://github.com/evolphotonics, August-October 2021) |
| 3 | +## Usage: sh Agreflmodelcalc.sh midul midgap midll CVul CVgap CVll [N] |
| 4 | + |
| 5 | +#!/bin/bash |
| 6 | + |
| 7 | +## Instead of hard coding, the scale thickness parameters for each scale type are passed as command line arguments (see usage) |
| 8 | + |
| 9 | +## upper lamina |
| 10 | +#ll=(56) ##FW Ag |
| 11 | +#ul=(64) ##FW Gland Ag |
| 12 | +#ul=(76) ##HW Gray |
| 13 | +#ul=(83) ##HW Ag |
| 14 | +#ul=(120) ##HW Coupling |
| 15 | +#midul="$(seq -s ' ' 0.001 5 501)" ## varying ul from 0 to 300 nm |
| 16 | +midul=$1 |
| 17 | + |
| 18 | +## air gap/lumen layer |
| 19 | +#gap=(1044) ##FW Ag |
| 20 | +#gap=(1420) ##FW Gland AG |
| 21 | +#gap=(819) ##HW Gray |
| 22 | +#gap=(1069) ##HW AG |
| 23 | +#gap=(1223) ##HW Coupling |
| 24 | +#midgap="$(seq -s ' ' 0.001 20 2001)" ## varying air gap from 0 to 2000 nm |
| 25 | +midgap=$2 |
| 26 | + |
| 27 | +## lower lamina |
| 28 | +#ll=(69) ##FW Ag |
| 29 | +#ll=(81) ##FW Gland Ag |
| 30 | +#ll=(83) ##HW Gray |
| 31 | +#ll=(67) ##HW Ag |
| 32 | +#ll=(105) ##HW Coupling |
| 33 | +#midll="$(seq -s ' ' 0.001 5 501)" ## varying ul from 0 to 500 nm |
| 34 | +midll=$3 |
| 35 | + |
| 36 | +## Coefficient of Variation |
| 37 | +#CV= #(0.08) (0.14) (0.07) ##FW Ag |
| 38 | +#CV= #(0.13) (0.10) (0.12) ##FW Gland Ag |
| 39 | +#CV= #(0.23) (0.14) (0.13) ##HW Gray |
| 40 | +#CV= #(0.11) (0.17) (0.14) ##HW Ag |
| 41 | +#CV= #(0.17) (0.28) (0.24) ##HW Coupling |
| 42 | +CVul=$4 |
| 43 | +CVgap=$5 |
| 44 | +CVll=$6 |
| 45 | + |
| 46 | +## Sample size |
| 47 | +if [ $7 -ge 0 ]; |
| 48 | +then |
| 49 | + N=$7 |
| 50 | +else |
| 51 | + N=200 |
| 52 | +fi |
| 53 | + |
| 54 | +count=1 ## counter start |
| 55 | +extension=txt ## for output file name |
| 56 | +suffix=e-9 ## for nanometer |
| 57 | + |
| 58 | +runfs () { |
| 59 | + local i=$1 |
| 60 | + local j=$2 |
| 61 | + local k=$3 |
| 62 | + |
| 63 | + |
| 64 | + ## increment and run Freesnell by passing the clparams as cl-arguments |
| 65 | + clparamul=$(echo $i$suffix) #"$i" |
| 66 | + clparamgap=$(echo $j$suffix) #"$j" |
| 67 | + clparamll=$(echo $k$suffix) #"$k" |
| 68 | + |
| 69 | + ## append direct output to dev null for speed |
| 70 | + /usr/local/bin/scm -v -f BanySilver.scm $clparamul $clparamgap $clparamll $i-$j-$k.$extension > /dev/null 2>&1 |
| 71 | + |
| 72 | + ## run verbosely by uncommenting the next line |
| 73 | + #echo "done with step $count"; |
| 74 | +} |
| 75 | + |
| 76 | +## We use the well known trick to average 3 uniform distributions to get a nearly normal distribution (See http://www.johndcook.com/blog/2009/02/12/sums-of-uniform-random-values/) |
| 77 | + |
| 78 | +## generate the normally distributed upper lamina values to average over |
| 79 | +#ul=$(for i in $(seq 1 $N) ; do awk -v seed=$RANDOM 'BEGIN{srand(seed); rdm=(((rand()*2)-1)+((rand()*2)-1)+((rand()*2)-1)); print (rdm*'"$CVul"'*'"$m"'+'"$m"')}'; done) |
| 80 | +if (( $(echo "$CVul > 0" |bc -l) )); ## Do only if the upper lamina layer thickness varies |
| 81 | +then |
| 82 | +uavg=0 ## initialize |
| 83 | +tol=5 ## tolerance of distribution mean from desired value, in nm |
| 84 | +while (( $(echo "$midul $uavg $tol" | awk '{print (sqrt(($1 - $2)^2) > $3)}') )) ## we shouldn't need this loop, but just in case |
| 85 | + do |
| 86 | + echo "$(for i in $(seq 1 $N) ; do awk -v seed=$RANDOM 'BEGIN{srand(seed); rdm=(((rand()*2)-1)+((rand()*2)-1)+((rand()*2)-1)); print (rdm*'"$CVul"'*'"$midul"'+'"$midul"')}'; done)" > ul |
| 87 | + uavg="$(awk '{sum+=$0} END { print sum/NR}' ul)" |
| 88 | + done; |
| 89 | +ul="$(cat ul)" |
| 90 | +else |
| 91 | +ul=$midul |
| 92 | +fi |
| 93 | + |
| 94 | +## generate the normally distributed air gap values to average over |
| 95 | +#gap=$(for i in $(seq 1 $N) ; do awk -v seed=$RANDOM 'BEGIN{srand(seed); rdm=(((rand()*2)-1)+((rand()*2)-1)+((rand()*2)-1)); print (rdm*'"$CVgap"'*'"$m"'+'"$m"')}'; done) |
| 96 | +if (( $(echo "$CVgap > 0" |bc -l) )); ## Do only if the air gap layer thickness varies |
| 97 | +then |
| 98 | +gavg=0 ## initialize |
| 99 | +tol=5 ## tolerance of distribution mean from desired value, in nm |
| 100 | +while (( $(echo "$midgap $gavg $tol" | awk '{print (sqrt(($1 - $2)^2) > $3)}') )) ## we shouldn't need this loop, but just in case |
| 101 | + do |
| 102 | + echo "$(for i in $(seq 1 $N) ; do awk -v seed=$RANDOM 'BEGIN{srand(seed); rdm=(((rand()*2)-1)+((rand()*2)-1)+((rand()*2)-1)); print (rdm*'"$CVgap"'*'"$midgap"'+'"$midgap"')}'; done)" > gap |
| 103 | + gavg="$(awk '{sum+=$0} END { print sum/NR}' gap)" |
| 104 | +done; |
| 105 | +gap="$(cat gap)" |
| 106 | +else |
| 107 | +gap=$midgap |
| 108 | +fi |
| 109 | + |
| 110 | +## generate the normally distributed lower lamina values to average over |
| 111 | +#ll=$(for i in $(seq 1 $N) ; do awk -v seed=$RANDOM 'BEGIN{srand(seed); rdm=(((rand()*2)-1)+((rand()*2)-1)+((rand()*2)-1)); print (rdm*'"$CVll"'*'"$m"'+'"$m"')}'; done) |
| 112 | +if (( $(echo "$CVll > 0" |bc -l) )); ## Do only if the lower lamina layer thickness varies |
| 113 | +then |
| 114 | +uavg=0 ## initialize |
| 115 | +tol=5 ## tolerance of distribution mean from desired value, in nm |
| 116 | +while (( $(echo "$midll $uavg $tol" | awk '{print (sqrt(($1 - $2)^2) > $3)}') )) ## we shouldn't need this loop, but just in case |
| 117 | + do |
| 118 | + echo "$(for i in $(seq 1 $N) ; do awk -v seed=$RANDOM 'BEGIN{srand(seed); rdm=(((rand()*2)-1)+((rand()*2)-1)+((rand()*2)-1)); print (rdm*'"$CVll"'*'"$midll"'+'"$midll"')}'; done)" > ll |
| 119 | + uavg="$(awk '{sum+=$0} END { print sum/NR}' ll)" |
| 120 | +done; |
| 121 | +ll="$(cat ll)" |
| 122 | +else |
| 123 | +ll=$midll |
| 124 | +fi |
| 125 | + |
| 126 | + |
| 127 | +## changing parameters simultaneously - if any of them are constant, then the corresponding for loop is only executed once |
| 128 | +for i in ${ul[@]}; do for j in ${gap[@]}; do for k in ${ll[@]}; do |
| 129 | + |
| 130 | + num_procs=100 ## number of concurrent jobs -- adjust according to your computer hardware specs |
| 131 | + num_jobs="\j" ## The prompt escape for number of jobs currently running |
| 132 | + while (( ${num_jobs@P} >= num_procs )); do |
| 133 | + wait ## wait before starting any more jobs |
| 134 | + done |
| 135 | + |
| 136 | + ## run each loop iteration in parallel |
| 137 | + runfs "$i" "$j" "$k" "$count" & |
| 138 | + |
| 139 | + ## increment counter |
| 140 | + count=$((count+1)) |
| 141 | + |
| 142 | +done; |
| 143 | +done; |
| 144 | +done; |
| 145 | + |
| 146 | + |
| 147 | +num_jobs="\j" ## The prompt escape for number of jobs currently running |
| 148 | +while (( ${num_jobs@P} > 0 )); do |
| 149 | + # echo "waiting for job" ${num_jobs@P} |
| 150 | + wait; ## for jobs to finish |
| 151 | +done |
| 152 | + |
| 153 | +## First calculate average spectra |
| 154 | +awk -v avgfile="avg.dat" -f avg.awk *.txt |
| 155 | + |
| 156 | +## Then calculate standard deviation of spectra |
| 157 | +awk -v avgfile="avg.dat" -f sd.awk *.txt |
| 158 | + |
| 159 | +## post processing/cleanup |
| 160 | +rm *.txt |
| 161 | +rm avg.dat |
| 162 | +mv avgwsd.dat $midul"_"$midgap"_"$midll"_"$CVul"_"$CVgap"_"$CVll"_"$N"_avg.dat" ## a more meaningful and unique name for the result |
| 163 | + |
| 164 | +exit 0 |
0 commit comments