Hi,
I find GMT pscoupe a bit tricky to use, and am thus trying to plot focal mechanisms and moment tensors from the side for profiles by first rotating them myself, and then using GMT psmeca.
On the left in the attached gif movie are example focal mechanisms with their r-theta-phi moment tensors next to them, plotted simply with psmeca in a Cartesian project. I understand those to be lower hemisphere projections, viewed from above as is the convention.
On the right are my attempts to rotate the same tensors into an x = along profile z = depth plane, as pscoupe can, and then trying to view them from the side, with a vertical (dip = 90 deg) plane.
The gif movie (use some viewer to stop and go back and forth) rotates this projection plane from the original N-S orientation (blue line) CW, and the red vector is meant to indicate how we are looking onto the plane. The right plot has the same events in the same locations as on the left, but after I rotate the matrix like I think it should be rotated.
What I find is that the psmeca projections on the side are jumping around, as if the routine was having a hard time with angles? Perhaps I am also plugging in things wrongly, but don’t know why/how.
I recall that psmeca had some issues before but those were fixed. Maybe some others remain?
I also attach the script used to make those plots.
Thanks much in advance!
#!/bin/bash
mode=2
if [ $mode = 1 ];then
head -15 $datadir/cmt/all.momt > tmp.dat
elif [ $mode = 2 ];then
cat <<EOF > tmp.dat
0 0 0 0 1 -1 0 0 0 23 0 0 1
0 0 0 0 -1 1 0 0 0 23 0 0 1
0 0 0 1 0 -1 0 0 0 23 0 0 1
0 0 0 -1 0 1 0 0 0 23 0 0 1
0 0 0 1 -1 0 0 0 0 23 0 0 1
0 0 0 -1 1 0 0 0 0 23 0 0 1
0 0 0 0 0 0 1 0 0 23 0 0 1
0 0 0 0 0 0 0 1 0 23 0 0 1
0 0 0 0 0 0 .99 .01 0 23 0 0 1
0 0 0 0 0 0 .01 .99 0 23 0 0 1
0 0 0 0 0 0 0 0 1 23 0 0 1
0 0 0 1 1 -2 0 0 0 23 0 0 1
0 0 0 -1 -1 2 0 0 0 23 0 0 1
EOF
head -2 $datadir/cmt/all.momt >> tmp.dat
fi
gawk 'BEGIN{y=5;x=0;}{x++;if(x==4){y--;x=1;}printf("%g %g ",x,y);
for(i=3;i<=NF;i++)printf("%s ",$i);printf("\n");}' tmp.dat > tmp.o.dat
cat tmp.o.dat
reg=-R0.5/3.75/0.5/5.5
proj=-Jx3
azi=0;gf=""
while [ $azi -le 360 ];do
echo $azi
cat tmp.o.dat | momt_project $azi > tmp.p.dat
#cat tmp.p.dat
((azim90 = azi - 90))
fs=""
for t in tmp.o tmp.p ;do
if [ $t = tmp ];then
l1=x;l2=y
else
l1=x;l2=z
fi
if [ $t = tmp.p2 ];then
gawk '{print($1,0,6-$2,$4,$5,$6,$7,$8,$9,$10,$11,$12,$13)}' tmp.o.dat | \
pscoupe -R0/4/0/5 -Aa0/0/90/10 -P -K -Sm1.5/-1 $proj > $t.ps
else
psmeca $proj $reg -Sm1.5/-1 $t.dat -Ba1:"$l1":/a1:"$l2": -K -P > $t.ps
fi
# 4 5 6 7 8 9
# X Y depth mrr mtt mff mrt mrf mtf
gawk 'BEGIN{xoff=0.3;yoff=0.1;dx=0.09;dy=0.075}{
m[1*3+1] = $4;m[1*3+2] =$7 ;m[1*3+3] = $8;
m[2*3+1] = $7;m[2*3+2] = $5;m[2*3+3] = $9;
m[3*3+1] = $8;m[3*3+2] =$9 ;m[3*3+3] = $6;
for(i=1;i<=3;i++)
for(j=1;j<=3;j++)
printf("%g %g 8 0 29 CM %4.1f\n",$1+xoff+(i-1)*dx,$2+yoff-(j-1)*dy,m[j*3+i]);
}' $t.dat | \
pstext $proj $reg -K -O >> $t.ps
if [ $t = tmp.o ];then
echo 1.5 4.5 $azi 1 | psxy $reg $proj -SV0.1/0.01/0.01 -Gblue -O -K >> $t.ps
echo 1.5 4.5 $azim90 1 | psxy $reg $proj -SV0.1/0.15/0.1 -Gred -O >> $t.ps
else
echo 100 100 0.1 | psxy $reg $proj -O -Sc.1 >> $t.ps
fi
modifybb $t.ps
#gv $t.ps ;exit
fs="$fs $t.ps"
done
epsmerge -par -x 2 -y 1 --orientation Landscape $fs > tmp.$azi.ps
#gv tmp.$azi.ps; exit
azil=`echo $azi | gawk '{printf("%010f",$1)}'`
/usr/bin/convert -rotate 90 -density 150 -background white -flatten -trim +repage tmp.$azi.ps tmp.$azil.gif &
gf="$gf tmp.$azil.gif"
rm tmp.p.dat
((azi=azi+10))
done
wait
gifsicle $gf > test_meca.gif
rm tmp*
