gfortran -o f90_fem_plnt f90_fem_plnt.f90

./f90_fem_plnt inp_4nod_ps100000.csv out_4nod_ps100000.csv
./f90_fem_plnt inp_4nod_ps050000.csv out_4nod_ps050000.csv
./f90_fem_plnt inp_4nod_ps020000.csv out_4nod_ps020000.csv
./f90_fem_plnt inp_4nod_ps010000.csv out_4nod_ps010000.csv
./f90_fem_plnt inp_4nod_ps005000.csv out_4nod_ps005000.csv
./f90_fem_plnt inp_4nod_ps002000.csv out_4nod_ps002000.csv
./f90_fem_plnt inp_4nod_ps001000.csv out_4nod_ps001000.csv
./f90_fem_plnt inp_4nod_ps000500.csv out_4nod_ps000500.csv
./f90_fem_plnt inp_4nod_ps000200.csv out_4nod_ps000200.csv
./f90_fem_plnt inp_4nod_ps000100.csv out_4nod_ps000100.csv
./f90_fem_plnt inp_4nod_ps000050.csv out_4nod_ps000050.csv
./f90_fem_plnt inp_4nod_ps000020.csv out_4nod_ps000020.csv
./f90_fem_plnt inp_4nod_ps000010.csv out_4nod_ps000010.csv
./f90_fem_plnt inp_4nod_ps000005.csv out_4nod_ps000005.csv
./f90_fem_plnt inp_4nod_ps000002.csv out_4nod_ps000002.csv
./f90_fem_plnt inp_4nod_ps000001.csv out_4nod_ps000001.csv

cat << EOT > _awk.awk
BEGIN{FS=",";s=0;r=0;E=ARGV[2];delete ARGV[2]}
NR==514,NR==533{s=s+\$8}
NR==336,NR==356{r=r+sqrt((\$2+\$4)**2+(\$3+\$5)**2)-sqrt((\$2)**2+(\$3)**2)}
END{s=s/20;d=r/21;print E,s,d}
EOT
awk -f _awk.awk out_4nod_ps000001.csv 1 > _sd.txt
awk -f _awk.awk out_4nod_ps000002.csv 2 >> _sd.txt
awk -f _awk.awk out_4nod_ps000005.csv 5 >> _sd.txt
awk -f _awk.awk out_4nod_ps000010.csv 10 >> _sd.txt
awk -f _awk.awk out_4nod_ps000020.csv 20 >> _sd.txt
awk -f _awk.awk out_4nod_ps000050.csv 50 >> _sd.txt
awk -f _awk.awk out_4nod_ps000100.csv 100 >> _sd.txt
awk -f _awk.awk out_4nod_ps000200.csv 200 >> _sd.txt
awk -f _awk.awk out_4nod_ps000500.csv 500 >> _sd.txt
awk -f _awk.awk out_4nod_ps001000.csv 1000 >> _sd.txt
awk -f _awk.awk out_4nod_ps002000.csv 2000 >> _sd.txt
awk -f _awk.awk out_4nod_ps005000.csv 5000 >> _sd.txt
awk -f _awk.awk out_4nod_ps010000.csv 10000 >> _sd.txt
awk -f _awk.awk out_4nod_ps020000.csv 20000 >> _sd.txt
awk -f _awk.awk out_4nod_ps050000.csv 50000 >> _sd.txt
awk -f _awk.awk out_4nod_ps100000.csv 100000 >> _sd.txt

# Theoritical values
cat << EOT > _awk.awk
BEGIN{
    Eg=ARGV[1]
    P=10
    Es=206000
    Ec=20600
    t=40
    rs=2300
    rc=3300
    nu=0.3
    lambda=1/(1+Es/Ec*t/rs*log(rc/rs)+Es/Eg*(1+nu)*t/rs)
    sigs=P*rs/t*(1-lambda)
    ds=P*(1-lambda)*rs*rs/Es/t
    print Eg,sigs,ds
}
EOT
awk -f _awk.awk 1 > _tsd.txt
awk -f _awk.awk 2 >> _tsd.txt
awk -f _awk.awk 5 >> _tsd.txt
awk -f _awk.awk 10 >> _tsd.txt
awk -f _awk.awk 20 >> _tsd.txt
awk -f _awk.awk 50 >> _tsd.txt
awk -f _awk.awk 100 >> _tsd.txt
awk -f _awk.awk 200 >> _tsd.txt
awk -f _awk.awk 500 >> _tsd.txt
awk -f _awk.awk 1000 >> _tsd.txt
awk -f _awk.awk 2000 >> _tsd.txt
awk -f _awk.awk 5000 >> _tsd.txt
awk -f _awk.awk 10000 >> _tsd.txt
awk -f _awk.awk 20000 >> _tsd.txt
awk -f _awk.awk 50000 >> _tsd.txt
awk -f _awk.awk 100000 >> _tsd.txt

cat << EOT > _inp.txt
10 $ 6 $ tr $ 0.3 $ 0.3
10 $ 0.5 $ 6
Theory     $ 1   $ ---   $ ---  $ ---
Calculated $ --- $ C0.2  $ 0    $ ---
EOT
awk -f /Users/kk/DATA_GMT/awk_klegend.awk _inp.txt > _legend.txt

range=1/100000/0/800
scale=10l/6
xga=a1g3
yga=g100a100
xlabel='Elastic Modulus of Bedrock Eg (MPa)'
ylabel='Stress of Penstock (MPa)'
fig=fig_str.eps
gmt set FONT_ANNOT_PRIMARY 12
gmt set MAP_ANNOT_OFFSET_PRIMARY 0.3c
gmt set FONT_LABEL 12p
gmt set MAP_LABEL_OFFSET 0.5c
gmt set MAP_TICK_LENGTH_PRIMARY 0.0c
gmt psbasemap -R$range -JX$scale -Bx$xga+l"$xlabel" -By$yga+l"$ylabel" -BWSen -P -K > $fig
awk '{print $1,$2}' _tsd.txt | gmt psxy -R -J -W1 -K -O >> $fig
awk '{print $1,$2}' _sd.txt | gmt psxy -R -J -SC0.2 -G#000000 -N -K -O >> $fig
source _legend.txt

range=1/100000/0/8
scale=10l/6
xga=a1g3
yga=g1a1
xlabel='Elastic Modulus of Bedrock Eg (MPa)'
ylabel='Displacement of Penstock (mm)'
fig=fig_dis.eps
gmt set FONT_ANNOT_PRIMARY 12
gmt set MAP_ANNOT_OFFSET_PRIMARY 0.3c
gmt set FONT_LABEL 12p
gmt set MAP_LABEL_OFFSET 0.5c
gmt set MAP_TICK_LENGTH_PRIMARY 0.0c
gmt psbasemap -R$range -JX$scale -Bx$xga+l"$xlabel" -By$yga+l"$ylabel" -BWSen -P -K > $fig
awk '{print $1,$3}' _tsd.txt | gmt psxy -R -J -W1 -K -O >> $fig
awk '{print $1,$3}' _sd.txt | gmt psxy -R -J -SC0.2 -G#000000 -N -K -O >> $fig
source _legend.txt

rm _*