fltdat2cquad

#!/bin/ksh
#
# reads flt patch geometry and slip from flt.dat and geom.in
# and writes to a geomview cquad file
#
# $Id: fltdat2cquad,v 1.1 2002/12/03 17:24:55 becker Exp becker $
#
# types of output: 
# 1: absolute slip value
# 2: normal stress
# 3: shear (strike) stress
type=${1-1} 
# cutoff, fraction of max
coff=${2-0.25}
# geometry
gfile=${3-geom.in}
# flt slip and stress data
fsfile=${4-flt.dat}
# fix max value?
fixmax=${5-0}

if [[ ! -s $gfile || ! -s $fsfile ]];then
    echo $0: $gfile or $fsfile not found
    exit
fi
tmpn=/tmp/$USER.$HOST.$$.pfd
trap "rm -f $tmpn.* ; exit" 0 1 2 15
#
# extract corners and vectors
#
cat $gfile | patch2xyzvec > $tmpn.cvec 2> /dev/null
# extract values: slip: strike, dip, normal absolute  
# stress: strike dip normal
gawk '{if((substr($1,1,1)!="#")&&(NF>8))\
    print($6,$8,$7,sqrt($6*$6+$8*$8+$7*$7),$9,$11,$10)}' \
    $fsfile > $tmpn.val
# determine extrema of values
minmax  -C $tmpn.val > $tmpn.ex
read usmin usmax udmin udmax unmin unmax uamin uamax ssmin ssmax sdmin sdmax snmin snmax < $tmpn.ex
# determine mean of values
gawk -f meanallcol.awk $tmpn.val > $tmpn.mean
read usmean udmean unmean uamean ssmean sdmean snmean < $tmpn.mean

# check file length
npatch=`lc $tmpn.cvec`
if [[ $npatch -ne `lc $tmpn.val` ]];then
    echo $0: length mismatch
    echo $0: cvec: `lc $tmpn.cvec`
    echo $0: slip: `lc $tmpn.val`
    exit
else
    echo $0: $npatch patches
fi
if [ $type -eq 1 ];then # absolute slip
    name=uabs
    use=25
    min=$uamin;max=$uamax;mean=$uamean
    label="u@-abs@-"
#    ls=a.5f.05
    ls=a1f.1
elif [ $type -eq 2 ];then # normal stress
    name=snorm
    use=28
    min=$snmin;max=$snmax;mean=$snmean
    label="@~s@~@-n@-"
    ls=a.2f.02
elif [ $type -eq 3 ];then # shear stress
    name=sstr
    use=26
    min=$ssmin;max=$ssmax;mean=$ssmean
    label="@~s@~@-s@-"
    ls=a.2f0.02
else
    echo $0: type $type not implemented
    exit
fi
#
# readjust min and max
#
echo $min $max $coff $mean | \
    gawk '{r=$2-$1;mi=$4-r*$3;ma=$4+r*$3;if(mi<$1)mi=$1;if(ma>$2)ma=$2;
	print(mi,ma)}' > \
    $tmpn.rem
read min max < $tmpn.rem
if [ $fixmax -ne 0 ];then
    min=0;max=$fixmax
fi

echo $0: type: $type min: $min max: $max
#
# output
#
paste $tmpn.cvec $tmpn.val | \
# cvec x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4 sx sy sz dx dy dz nx ny nz  ustrike udip unormal uabs sstrike sdip snormal
#      1        4        7        10       13       16       19    21    22     23    24     25    26     27    28
gawk --assign use=$use --assign min=$min --assign max=$max \
'BEGIN{
 range=max-min;\
 print("CQUAD");\
}\
{\
 x=($(use)-min)/range;\
 col[1]=x;col[2]=1-((x-.5)**2)*4.;col[3]=1.0-x;\
 for(i=1;i<=3;i++){\
   if(col[i]<0)col[i]=0.0;\
   if(col[i]>1)col[i]=1.0;\
 }\
 for(i=1;i<=4;i++)
    printf("%g %g %g  %g %g %g 1.0\n",
	$((i-1)*3+1),$((i-1)*3+2),$((i-1)*3+3),
	col[1],col[2],col[3]);
}' > flt.$name.cquad
#
# colorbar
#
dy=0.05
gawk --assign min=$min --assign max=$max --assign steps=150 \
    --assign dy=$dy \
'BEGIN{\
 range=max-min;\
 dx=range/steps;\
 xs=min+0.5*dx;
 for(x=min;x<max;x+=dx){\
  xp=(xs-min)/range;\
  col[1]=xp;col[2]=1-((xp-.5)**2)*4.;col[3]=1.0-xp;\
  for(i=1;i<=3;i++){\
   col[i] *= 255;
   if(col[i]<0)col[i]=0;\
   if(col[i]>255)col[i]=255;\
  }\
  printf("> -G%i/%i/%i \n",col[1],col[2],col[3]);\
  print(x,0);print(x+dx,0);print(x+dx,dy);print(x,dy);
  xs+=dx;
 }\
}' > $tmpn.cb
psxy $tmpn.cb -R$min/$max/0.001/$dy -JX4.6/.5 \
    -B$ls/f100:.$label:weSn -P -N -M  > flt.$name.cb.ps
modifybb flt.$name.cb.ps 65 40 405 170 > /dev/null
bbtofront flt.$name.cb.ps 
echo $0: output in flt.$name.cquad and flt.$name.cb.ps