#!/bin/bash
# 
# test different ways of extracting the solutions
#
layers=${1--1}			# -1: all layers
runcode=${2-1}			# run the code
cspec=${3-1}			# compare spectral, rather than spatial solutions
csol=${3-0}			# compare solution extraction methods for new code

lmax=31				# max degree of spectral comparison


model=2
if [ $model -eq 1 ];then
# A, plate motion only 
    old_model=$datadir/flow_field/finite_strain/pmAd/;dens_scale=0.0;free_slip="";visc=A;lmax=31
elif [ $model -eq 2 ];then
# A, free slip density driven
    old_model=$datadir/flow_field/finite_strain/pmAsmean_nt_free/;dens_scale=0.2;free_slip="-fs";visc=A;lmax=31
elif [ $model -eq 3 ];then
# F, plate motion only
    old_model=$datadir/flow_field/finite_strain/pmFd/;dens_scale=0.0;free_slip="";visc=F;lmax=63
elif [ $model -eq 4 ];then
# F, free slip
    old_model=$datadir/flow_field/finite_strain/pmFsmean_nt_free/;dens_scale=0.2;free_slip="-fs";visc=F;lmax=31
elif [ $model -eq 5 ];then
# F, plate motion and density
    old_model=$datadir/flow_field/finite_strain/pmFsmean_nt/;dens_scale=0.2;free_slip="";visc=F
fi


reg=-R0/360/-90/90


if [ `echo $layers | gawk '{if((NF==1)&&($1==-1))print(1);else print(0)}'` -eq 1 ];then
    # all layers 
    n=`lc $old_model/vdepth.dat`
    layers=`gawk --assign n=$n 'BEGIN{for(i=1;i<=n;i++)printf("%i ",i)}'`

fi
echo $0: checking layers $layers

rm  lrlog.?.$model.dat 2> /dev/null


if [ $runcode -eq 1 ];then
#
# prepare input files
#
    # convert plate velocities
    cat $HOME/plates/rick_expansions/my.nuvel.$lmax.1.0.cs | abconvert 7 0 1 6 | \
	gawk '{if(NR==1)print($1,0,0,1,2,0);else print($0)}' > pvel.sh.dat
    # convert density models (this would be OK, but 31...63 issue)
    convert_she_model $datadir/tomography/smean_nt/d.31.m.ab  > dens.sh.dat
    rm psol.dat tsol.dat vsol.*.bin 2> /dev/null
    # run code
    hc -pptsol -ds $dens_scale $free_slip -vf visc.$visc
fi
for layer in $layers;do

    if [ $cspec -eq 1 ];then

#
    # spectral comparison
#
# pol/tor
	hc_extract_sh_layer sol.bin $layer | gawk '{if(NR==1)print($1);else print($1,$2)}' > sh.orig.p 2> /dev/null
	hc_extract_sh_layer sol.bin $layer | gawk '{if(NR==1)print($1);else print($3,$4)}' > sh.orig.t 2> /dev/null
# radial
	hc_extract_sh_layer sol.bin $layer 1 | gawk '{if(NR==1)print($1);else print($0)}' > sh.orig.r 2> /dev/null

	for t in vp vt vr;do
	    if [ ! -s $old_model/$t.$layer.grd ];then
		echo $0: $old_model/$t.$layer.grd not found
		exit
	    fi
	done
	rm sh.old.? 2> /dev/null
	shana $lmax  $old_model/vr.$layer > sh.old.r  2> /dev/null
	cp $old_model/vp.$layer.grd vec_p.grd
	cp $old_model/vt.$layer.grd vec_t.grd
	shana $lmax vec_p > sh.old.pt 2> /dev/null
	rm vec_p.grd vec_t.grd
	gawk '{if(NR==1)print($1);else print($1,$2)}'  sh.old.pt > sh.old.p
	gawk '{if(NR==1)print($1);else print($3,$4)}'  sh.old.pt > sh.old.t
	rm sh.old.pt

    # compute correlation and linear regression
	for t in r p t;do
	    if [[ ! -s sh.orig.$t || ! -s sh.old.$t ]];then
		echo $0: sh.orig.$t or sh.old.$t not found
		exit
	    fi
	    cat sh.orig.$t  sh.old.$t | abconvert 5 | head -5
	    echo  `oneline $layer vdepth.dat` `paste sh.old.$t sh.orig.$t | gawk '{if(NR>1)print($1,$3);print($2,$4)}' | gawk -f linreg.awk` `cat sh.orig.$t  sh.old.$t | abconvert 6` >> lrlog.$t.$model.dat
	    tail -1 lrlog.$t.$model.dat
	done
    else
#
    # spatial
#
	if [ $csol -eq 1 ];then
#
# extract the solution as spherical harmonics directly from file
#
# pol/tor
	    hc_extract_sh_layer sol.bin $layer > sh.orig.pt 2> /dev/null
# radial
	    hc_extract_sh_layer sol.bin $layer 1 > sh.orig.r 2> /dev/null
#
#
# convert to spatial using HC type routines
#
	    cat sh.orig.pt | sh_syn > rspectral.vtvp 2> /dev/null
	    cat sh.orig.r  | sh_syn > rspectral.vr 2> /dev/null
# grid
	    
	    surface $reg -I2 rspectral.vr -Grspectral.vr.grd
	    gawk '{print($1,$2,$3)}' rspectral.vtvp | \
		surface $reg -I2 -Grspectral.vt.grd
	    gawk '{print($1,$2,$4)}' rspectral.vtvp | \
		surface $reg -I2 -Grspectral.vp.grd
#
# reexpand using HC routines
	    cat rspectral.vr   | sh_ana 63 0 >  sh.reexpanded.r 2> /dev/null
	    cat rspectral.vtvp | sh_ana 63 1 >  sh.reexpanded.pt 2> /dev/null
	    
	    
	    rm rspectral.vtvp rspectral.vr 
	    
#
# reexpand using my old routines
#
	    shana 63 rspectral.vr > sh.mreexpanded.r 2> /dev/null
	    cp rspectral.vt.grd vec_t.grd
	    cp rspectral.vp.grd vec_p.grd
	    shana 63 vec_p > sh.mreexpanded.pt 2> /dev/null
	    rm vec_p.grd vec_t.grd
# extract using my old routines
	    gawk '{if(NR==1)print($1);else print($0)}' sh.orig.r | \
		shsyn 2 0 mspectral.vr 2> /dev/null
	    gawk '{if(NR==1)print($1);else print($0)}'  sh.orig.pt  | \
		shsyn 2 0 3 2> /dev/null
	    mv vec_p.grd mspectral.vp.grd
	    mv vec_t.grd mspectral.vt.grd
# plot 
	    for m in rspectral mspectral ;do
		for t in vr vt vp;do
		    grd2map2 $m.$t vel $t > /dev/null 2> /dev/null
		done
	    done
	    rm rspectral.v?.grd mspectral.v?.grd
	    
#
# comparison of different ways of extracting the solution
#
	    epsmerge -par --print --postscript --orientation Landscape -x 3 -y 3 \
		ospatial.vr.ps ospatial.vp.ps ospatial.vt.ps \
		rspectral.vr.ps rspectral.vp.ps rspectral.vt.ps \
		mspectral.vr.ps mspectral.vp.ps mspectral.vt.ps \
		> sol.vec.ps 2> /dev/null
	fi

#
# convert velocities from spatial solution output
#
	i=1
	for f in lon lat vr vt vp;do
	    if [ ! -s vsol.$layer.bin ];then
		echo $0: vsol.$layer.bin not found
		exit
	    fi
	    cat vsol.$layer.bin | extract_bin 5 $i 0 1 > tmp.$f 2> /dev/null
	    ((i=i+1))
	done
	for f in vr vt vp;do
	    paste tmp.lon tmp.lat tmp.$f | \
		surface $reg -I2 -Gospatial.$f.grd
	    rm tmp.$f
	    grd2map2 ospatial.$f vel $f > /dev/null 2> /dev/null
	done
	
#
# quantitative comparison
#
	for t in vr vt vp;do
	    if [ ! -s $old_model/$t.$layer.grd ];then
		echo $0: $old_model/$t.$layer.grd not found
		exit
	    fi
	    grd2map2 $old_model/$t.$layer vel $t > /dev/null 2> /dev/null
	    mv $old_model/$t.$layer.ps old.$t.ps
	    echo $0: solution compare: $t: `calc_grd_correlation $old_model/$t.$layer.grd ospatial.$t.grd -1 1`
	done
#
# comparison between old and new HC code
#
	epsmerge -par --print --postscript --orientation Landscape -x 3 -y 2 \
	    old.vr.ps old.vp.ps old.vt.ps \
	    ospatial.vr.ps ospatial.vp.ps ospatial.vt.ps \
	    > sol.compare.ps 2> /dev/null
	
	rm  old.vr.ps old.vp.ps old.vt.ps \
	    ospatial.vr.ps ospatial.vp.ps ospatial.vt.ps \
	    rspectral.vr.ps rspectral.vp.ps rspectral.vt.ps \
	    mspectral.vr.ps mspectral.vp.ps mspectral.vt.ps 2> /dev/null
    fi				# end spatial

done