[cig-commits] r16292 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Sat Feb 20 22:09:05 PST 2010
Author: becker
Date: 2010-02-20 22:09:05 -0800 (Sat, 20 Feb 2010)
New Revision: 16292
Modified:
mc/1D/hc/trunk/calc_vel_and_plot
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/prem_util.c
mc/1D/hc/trunk/sh_exp.c
Log:
Slight robustness improvements for sh_syn
Added proper velocity grid extraction for the calc_vel_and_plot example
script.
Modified: mc/1D/hc/trunk/calc_vel_and_plot
===================================================================
--- mc/1D/hc/trunk/calc_vel_and_plot 2010-02-21 05:18:15 UTC (rev 16291)
+++ mc/1D/hc/trunk/calc_vel_and_plot 2010-02-21 06:09:05 UTC (rev 16292)
@@ -49,7 +49,7 @@
for plot in $plots;do
- if [[ $plot -eq 1 || $plot -eq 2 ]];then # velocities or geoid
+ if [[ $plot -eq 1 || $plot -eq 2 || $plot -eq 4 ]];then # velocities or geoid
f=vel
hc_opt=""
@@ -62,154 +62,185 @@
if [ $compute -eq 1 ];then
rm $f.sol.bin 2> /dev/null
+ echo $0: recomputing
hc $verb -dens $dm $dt $hc_opt -ds $dfac $tbc -vf $vf 2> hc.log
if [ ! -s $f.sol.bin ];then
echo $0: error, no output $f.sol.bin produced
exit
fi
cat hc.log
+ else
+ echo $0: reusing $f.sol.bin
fi
+ if [ ! -s $f.sol.bin ];then
+ echo $0: error: $f.sol.bin not found
+ exit
+ fi
+
+ if [ $plot -le 3 ];then # actually plotting
+
#
# take care of GMT crap
#
- if [ -s $HOME/.gmtdefaults ];then
- gf=$HOME/.gmtdefaults
- elif [ -s $HOME/.gmtdefaults4 ];then
- gf=$HOME/.gmtdefaults4
- elif [ -s .gmtdefaults4 ];then
- gf=.gmtdefaults4
- elif [ -s .gmtdefaults ];then
- gf=.gmtdefaults
- fi
- if [ -s $gf ];then
- old_media=`grep PAPER_MEDIA $gf | gawk '{print($3)}'`
- fi
-
- gmtset PAPER_MEDIA letter+
+ if [ -s $HOME/.gmtdefaults ];then
+ gf=$HOME/.gmtdefaults
+ elif [ -s $HOME/.gmtdefaults4 ];then
+ gf=$HOME/.gmtdefaults4
+ elif [ -s .gmtdefaults4 ];then
+ gf=.gmtdefaults4
+ elif [ -s .gmtdefaults ];then
+ gf=.gmtdefaults
+ fi
+ if [ -s $gf ];then
+ old_media=`grep PAPER_MEDIA $gf | gawk '{print($3)}'`
+ fi
+
+ gmtset PAPER_MEDIA letter+
#
# plotting setting
#
- inc=1;
- greg=-R0/359/-89.5/89.5 # grid range
- proj=-JH180/7i;preg=-R0/360/-90/90 # plotting range
-
- if [ $plot -eq 1 ];then # plot the geoid
+ inc=1;
+ greg=-R0/359/-89.5/89.5 # grid range
+ proj=-JH180/7i;preg=-R0/360/-90/90 # plotting range
- ofile=geoid.ps # output PS filename
-
- echo $0: plotting geoid to $ofile
-
-
- makecpt -T-150/150/10 -Chaxby > $tmpn.cpt # makes colormap
-
+ if [ $plot -eq 1 ];then # plot the geoid
+
+ ofile=geoid.ps # output PS filename
+
+ echo $0: plotting geoid to $ofile
+
+
+ makecpt -T-150/150/10 -Chaxby > $tmpn.cpt # makes colormap
+
# go from spherical harmonics to spatial
- cat geoid.ab | sh_syn 0 0 $inc 2> /dev/null | \
- xyz2grd $greg -I$inc -G$tmpn.grd # and convert to netcdf GRD
-
+ cat geoid.ab | sh_syn 0 0 $inc 2> /dev/null | \
+ xyz2grd $greg -I$inc -G$tmpn.grd # and convert to netcdf GRD
+
# make postscript plot from grd file
- grdimage $tmpn.grd $preg $proj -P -C$tmpn.cpt -K > $ofile # first: -K
+ grdimage $tmpn.grd $preg $proj -P -C$tmpn.cpt -K > $ofile # first: -K
# coastline
- pscoast $proj $preg -A70000 -K -O -Dc -W2 >> $ofile # all intermediat: -O -K
+ pscoast $proj $preg -A70000 -K -O -Dc -W2 >> $ofile # all intermediat: -O -K
# colorbar
- psscale -N50 -D3.5i/-.3i/3i/.25ih -E -C$tmpn.cpt -B50/:"[m]": -O >> $ofile # end: -O
-
+ psscale -N50 -D3.5i/-.3i/3i/.25ih -E -C$tmpn.cpt -B50/:"[m]": -O >> $ofile # end: -O
+
# ps file done
-
- elif [[ $plot -eq 2 || $plot -eq 3 ]];then # plot velocities or tractions
-
+
+ elif [[ $plot -eq 2 || $plot -eq 3 ]];then # plot velocities or tractions
+
# extract all depths , use second colume, print depth in [km]
- hc_extract_sh_layer $f.sol.bin 2 4 | gawk '{print($2)}' > $dfile.dat
-
+ hc_extract_sh_layer $f.sol.bin 2 4 | gawk '{print($2)}' > $dfile.dat
+
# count number of lines (assign command output to variable)
- nl=`wc -l $dfile.dat | gawk '{print($1)}'`
-
- if [ $plot -eq 2 ];then
+ nl=`wc -l $dfile.dat | gawk '{print($1)}'`
+
+ if [ $plot -eq 2 ];then
# colormaps for radial flow
- makecpt -T-1.5/1.5/.125 -Chaxby > $tmpn.cpt
- scsp=.5 # scale spacing
- type="v";units="cm/yr"
- else # tractions
- scsp=50
- type="t";units="MPa"
- makecpt -T-150/150/10 -Chaxby > $tmpn.cpt
- fi
-
+ makecpt -T-1.5/1.5/.125 -Chaxby > $tmpn.cpt
+ scsp=.5 # scale spacing
+ type="v";units="cm/yr"
+ else # tractions
+ scsp=50
+ type="t";units="MPa"
+ makecpt -T-150/150/10 -Chaxby > $tmpn.cpt
+ fi
+
# vector spacing
- vinc=10;vreg=-R0/350/-85/85
-
+ vinc=10;vreg=-R0/350/-85/85
+
#
# loop through some layers
#
- lc=1
- for layer in $layers;do # loop
-
- if [ $layer -gt $nl ];then # check
- echo $0: error, layer $layer out of bounds $nl
- exit
- fi
-
+ lc=1
+ for layer in $layers;do # loop
+
+ if [ $layer -gt $nl ];then # check
+ echo $0: error, layer $layer out of bounds $nl
+ exit
+ fi
+
#
# extract depth
- z=`gawk '{if(NR==l)printf("%.0f",$1)}' l=$layer $dfile.dat `
-
- ofile=$f.$layer.ps # output file name
-
- echo $0: plotting $f at $z to $ofile
+ z=`gawk '{if(NR==l)printf("%.0f",$1)}' l=$layer $dfile.dat `
+
+ ofile=$f.$layer.ps # output file name
+
+ echo $0: plotting $f at $z to $ofile
#
# extract radial velocity of this layer
#
# extract layer $layer as spherical harmonics , radial component
# convert to spatial
# make netcdf GRD
- hc_extract_sh_layer $f.sol.bin $layer 1 0 | sh_syn 0 0 $inc 2> /dev/null | \
- xyz2grd $greg -I$inc -G$tmpn.grd
+ hc_extract_sh_layer $f.sol.bin $layer 1 0 | sh_syn 0 0 $inc 2> /dev/null | \
+ xyz2grd $greg -I$inc -G$tmpn.grd
#
# extract vx and vy velocities
#
- hc_extract_sh_layer $f.sol.bin $layer 2 0 | sh_syn 0 0 $vinc 2> /dev/null > $tmpn.dat
+ hc_extract_sh_layer $f.sol.bin $layer 2 0 | sh_syn 0 0 $vinc 2> /dev/null > $tmpn.dat
#
- gawk '{print($1,$2,$4)}' $tmpn.dat | xyz2grd -G$tmpn.vx.grd $vreg -I$vinc # vx = vphi, tphi
- gawk '{print($1,$2,-$3)}' $tmpn.dat | xyz2grd -G$tmpn.vy.grd $vreg -I$vinc # vy = -vtheta,-ttheta
-
+ gawk '{print($1,$2,$4)}' $tmpn.dat | xyz2grd -G$tmpn.vx.grd $vreg -I$vinc # vx = vphi, tphi
+ gawk '{print($1,$2,-$3)}' $tmpn.dat | xyz2grd -G$tmpn.vy.grd $vreg -I$vinc # vy = -vtheta,-ttheta
+
# computing vector length
- grdmath $tmpn.vx.grd $tmpn.vy.grd R2 SQRT = $tmpn.abs.grd
-
- if [ $lc -eq 1 ];then # first time around
+ grdmath $tmpn.vx.grd $tmpn.vy.grd R2 SQRT = $tmpn.abs.grd
+
+ if [ $lc -eq 1 ];then # first time around
#
# get scale from multiples of mean of first layer
- hmean=`grdinfo -C -L2 $tmpn.abs.grd | gawk '{printf("%.1f",$12)}'` # compute mean
- scale=`echo $hmean | gawk '{print($1*5)}'` # scale
- fi
+ hmean=`grdinfo -C -L2 $tmpn.abs.grd | gawk '{printf("%.1f",$12)}'` # compute mean
+ scale=`echo $hmean | gawk '{print($1*5)}'` # scale
+ fi
#
# actually plot velocities
#
- grdimage $tmpn.grd $preg $proj -P -C$tmpn.cpt -K > $ofile # make a colormap
+ grdimage $tmpn.grd $preg $proj -P -C$tmpn.cpt -K > $ofile # make a colormap
# coastline
- pscoast $proj $preg -A70000 -K -O -Dc -W2 >> $ofile
-
+ pscoast $proj $preg -A70000 -K -O -Dc -W2 >> $ofile
+
#
- echo 0.05 0 14 0 0 ML "z = $z km" | pstext -O -K -N -R0/1/0/1 -JX7/3.5i >> $ofile
- echo 0.75 0 14 0 0 ML "@~\341@~$type at -h@-@~\361@~ = $hmean $units" | \
- pstext -O -K -N -R0/1/0/1 -JX7/3.5i >> $ofile
+ echo 0.05 0 14 0 0 ML "z = $z km" | pstext -O -K -N -R0/1/0/1 -JX7/3.5i >> $ofile
+ echo 0.75 0 14 0 0 ML "@~\341@~$type at -h@-@~\361@~ = $hmean $units" | \
+ pstext -O -K -N -R0/1/0/1 -JX7/3.5i >> $ofile
# plot a vector field
- grdvector $tmpn.vx.grd $tmpn.vy.grd \
- -T $reg $proj -Q0.015i/0.12i/0.045in.2 -S$scale -O -K -G128/0/0 -W0.5 >> $ofile
+ grdvector $tmpn.vx.grd $tmpn.vy.grd \
+ -T $reg $proj -Q0.015i/0.12i/0.045in.2 -S$scale -O -K -G128/0/0 -W0.5 >> $ofile
# scale
- psscale -D50 -D3.5i/-.3i/3i/.25ih -E -C$tmpn.cpt \
- -B$scsp/:"$type at -r@- [$units]": -O >> $ofile
+ psscale -D50 -D3.5i/-.3i/3i/.25ih -E -C$tmpn.cpt \
+ -B$scsp/:"$type at -r@- [$units]": -O >> $ofile
+
+ ((lc=lc+1))
+ done
+ fi
+ if [ -s $gf ];then # resset GMT
+ gmtset PAPER_MEDIA $old_media
+
+ fi
+ else # only extracting grids
+ vinc=1 # spacing
+ vreg=-R0/360/-89.5/89.5
+ echo $0: writing layer depth to vdepth.dat
+ hc_extract_sh_layer $f.sol.bin 2 4 | gawk '{print($2)}' > vdepth.dat
+ nl=`wc -l $dfile.dat | gawk '{print($1)}'`
+ i=1
+ while [ $i -le $nl ];do
+ hc_extract_sh_layer $f.sol.bin $i 1 0 | sh_syn 0 0 0 360 -89.5 89.5 $vinc $vinc 2> /dev/null | \
+ xyz2grd $vreg -I$vinc -Gvr.$i.grd # radial component
+ # theta, phi
+ hc_extract_sh_layer $f.sol.bin $i 2 0 | sh_syn 0 0 0 360 -89.5 89.5 $vinc $vinc 2> /dev/null > $tmpn.dat
+ #
+ gawk '{print($1,$2,$3)}' $tmpn.dat | xyz2grd -Gvt.$i.grd $vreg -I$vinc # vx = vphi, tphi
+ gawk '{print($1,$2,$4)}' $tmpn.dat | xyz2grd -Gvp.$i.grd $vreg -I$vinc # vy = -vtheta,-ttheta
+ echo $0: extracting layer $i out of $nl
+ exit
- ((lc=lc+1))
+ ((i=i+1))
done
+
fi
- if [ -s $gf ];then # resset GMT
- gmtset PAPER_MEDIA $old_media
- fi
-
-
done
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2010-02-21 05:18:15 UTC (rev 16291)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2010-02-21 06:09:05 UTC (rev 16292)
@@ -189,7 +189,8 @@
undefined and FALSE else
*/
ggrd_boolean ggrd_grdtrack_interpolate_rtp(double r,double t,double p,
- struct ggrd_gt *g,double *value,
+ struct ggrd_gt *g,
+ double *value,
ggrd_boolean verbose,
ggrd_boolean shift_to_pos_lon)
{
@@ -218,6 +219,7 @@
x[2] = (1.0-r) * 6371.0; /* depth in [km] */
if(g->zlevels_are_negative) /* adjust for depth */
x[2] = -x[2];
+
result = ggrd_grdtrack_interpolate(x,TRUE,g->grd,g->f,
g->edgeinfo,g->mm,g->z,
g->nz,value,verbose,
@@ -777,25 +779,25 @@
*/
#ifndef USE_GMT3
ggrd_boolean ggrd_grdtrack_interpolate(double *in, /* lon/lat/z [2/3] in degrees/km */
- ggrd_boolean three_d, /* use 3-D inetrpolation or 2-D? */
- struct GRD_HEADER *grd, /* grd information */
- float *f, /* data array */
- struct GMT_EDGEINFO *edgeinfo, /* edge information */
- int mm, /* nx * ny */
- float *z, /* depth layers */
- int nz, /* number of depth layers */
+ ggrd_boolean three_d, /* use 3-D inetrpolation or 2-D? */
+ struct GRD_HEADER *grd, /* grd information */
+ float *f, /* data array */
+ struct GMT_EDGEINFO *edgeinfo, /* edge information */
+ int mm, /* nx * ny */
+ float *z, /* depth layers */
+ int nz, /* number of depth layers */
double *value, /* output value */
- ggrd_boolean verbose,
- struct GMT_BCR *loc_bcr)
+ ggrd_boolean verbose,
+ struct GMT_BCR *loc_bcr)
#else
ggrd_boolean ggrd_grdtrack_interpolate(double *in, /* lon/lat/z [2/3] in degrees/km */
- ggrd_boolean three_d, /* use 3-D inetrpolation or 2-D? */
- struct GRD_HEADER *grd, /* grd information */
- float *f, /* data array */
- struct GMT_EDGEINFO *edgeinfo, /* edge information */
- int mm, /* nx * ny */
- float *z, /* depth layers */
- int nz, /* number of depth layers */
+ ggrd_boolean three_d, /* use 3-D inetrpolation or 2-D? */
+ struct GRD_HEADER *grd, /* grd information */
+ float *f, /* data array */
+ struct GMT_EDGEINFO *edgeinfo, /* edge information */
+ int mm, /* nx * ny */
+ float *z, /* depth layers */
+ int nz, /* number of depth layers */
double *value, /* output value */
ggrd_boolean verbose,
struct BCR *loc_bcr
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2010-02-21 05:18:15 UTC (rev 16291)
+++ mc/1D/hc/trunk/hc.h 2010-02-21 06:09:05 UTC (rev 16292)
@@ -61,7 +61,9 @@
#define HC_PI M_PI
+#define HC_TWOPI 6.2831853071795864769252867665590
+
#define HC_SCALAR 0
#define HC_VECTOR 1
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2010-02-21 05:18:15 UTC (rev 16291)
+++ mc/1D/hc/trunk/hc_init.c 2010-02-21 06:09:05 UTC (rev 16292)
@@ -851,8 +851,10 @@
hc_vecrealloc(&hc->dvisc,hc->nradp2,"hc_assign_density");
hc->r[0] = hc->r_cmb; /* CMB */
- if(hc->rden[0] <= hc->r[0])
- HC_ERROR("hc_assign_density","first density layer has to be above internal CMD limit");
+ if(hc->rden[0] <= hc->r[0]){
+ fprintf(stderr,"hc_assign_density: rden[0]: %g r[0]: %g\n",hc->rden[0], hc->r[0]);
+ HC_ERROR("hc_assign_density","first density layer has to be above internal CMB limit");
+ }
for(i=0;i<hc->nrad;i++) /* density layers */
hc->r[i+1] = hc->rden[i];
if(hc->rden[hc->nrad-1] >= 1.0)
Modified: mc/1D/hc/trunk/prem_util.c
===================================================================
--- mc/1D/hc/trunk/prem_util.c 2010-02-21 05:18:15 UTC (rev 16291)
+++ mc/1D/hc/trunk/prem_util.c 2010-02-21 06:09:05 UTC (rev 16292)
@@ -79,7 +79,7 @@
*/
double prem_compute_dpval(double *y, double *par,
- int np, double fac)
+ int np, double fac)
{
int i;
double val;
@@ -106,7 +106,7 @@
*/
double prem_vs_voigt(double vsh, double vsv, double vph, double vpv,
- double eta)
+ double eta)
{
double v;
v = sqrt((1./15.) * ( (1.0-2.0*eta) * vph*vph + vpv*vpv +
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2010-02-21 05:18:15 UTC (rev 16291)
+++ mc/1D/hc/trunk/sh_exp.c 2010-02-21 06:09:05 UTC (rev 16292)
@@ -1413,10 +1413,10 @@
/* print coordinates */
if(!use_3d){
/* print lon lat */
- fprintf(out,"%14.7f %14.7f\t",lon,lat);
+ fprintf(out,"%12.3f %12.3f\t",lon,lat);
}else{
/* print lon lat z[i] */
- fprintf(out,"%14.7f %14.7f %14.7f\t",lon,lat,z);
+ fprintf(out,"%12.3f %12.3f %14.7f\t",lon,lat,z);
}
for(k=0;k<shps;k++) /* loop through all scalars */
fprintf(out,"%14.7e ",data[l+npoints*k]);
More information about the CIG-COMMITS
mailing list