# GMT求救
```
REM Original batch file was created with GMT4 by Wan-Chun Chung (鍾琬君) ,Department of Geosciences, National Taiwan University
REM Modified by Jyr-Ching Hu with GMT 5
REM Modified by Bo-Han Tsao and Jyr-Ching Hu with GMT 6 modules
REM gmt configuration
gmt gmtset FONT_TITLE=20p,AvantGarde-Demi,black
gmt gmtset FONT_LABEL = 12p,Helvetica,black
gmt gmtset FONT_ANNOT_PRIMARY = 9p,Helvetica,black
REM set variables
set prefix=EX06
set grd=C:\gmtfile\Lecture06\grdfiles\Taiwan40m_WGS84.nc
set outgrd=tw_cent.nc
set cpt1=topo_gray.cpt
set cpt2=fpcolor.cpt
set int=tw_cent.grd.int
set x1=119.5
set x2=122.8
set y1=21.4098756634
set y2=25.5000722922
set range=%x1%/%x2%/%y1%/%y2%
::process data
gmt grdcut %grd% -G%outgrd% -R%range% -V
gmt makecpt -Cgray -T0/4000/500 -I -N > topo_gray.cpt -V
gmt project -C120.40/24.10 -E121.45/23.65 -G.5 -Q > track.A
gmt grdtrack track.A -G%outgrd% | gawk "{ print $3, $4 }" > trackA.xyp
echo 100 0 >> trackA.xyp
echo 0 0 >> trackA.xyp
::plot
gmt begin %prefix% jpg A+m0.5c
gmt grdimage %outgrd% -R%range% -Jm121.0/1.5i -C%cpt1% -I+a60+ne1.0 -V
gmt coast -W1p,0 -S255/255/255 -Lg119.80/24.20+c15+w50+f+lkm -Bx1f0.5 -By1f0.5 -BWSne -Df -V
REM ****EQ classify with focal depth****
gawk "{print $1, $2, -$3, $4}" CWAEQs_1995-2022.txt > CWAEQs_1995-2022.gmt
gawk "{if ($3 >= 0 && $3 < 15) print $1, $2, $3, $4 }" CWAEQs_1995-2022.gmt > hypo_0_15.gmt
gawk "{if ($3 >= 15 && $3 < 30) print $1, $2, $3, $4 }" CWAEQs_1995-2022.gmt > hypo_15_30.gmt
gawk "{if ($3 >= 30 && $3 < 50) print $1, $2, $3, $4 }" CWAEQs_1995-2022.gmt > hypo_30_50.gmt
gawk "{if ($3 >= 50 ) print $1, $2, $3, $4 }" CWAEQs_1995-2022.gmt > hypo_50+.gmt
gawk "{print $1, $2, $4*0.01}" hypo_0_15.gmt | gmt plot -Sc -W0.00008i,255/106/106 -V
gawk "{print $1, $2, $4*0.01}" hypo_15_30.gmt | gmt plot -Sc -W0.00008i,250/250/0 -V
gawk "{print $1, $2, $4*0.01}" hypo_30_50.gmt | gmt plot -Sc -W0.00008i,50/250/100 -V
gawk "{print $1, $2, $4*0.01}" hypo_50+.gmt | gmt plot -Sc -W0.00008i,0/191/255 -V
REM ***** Basement high *****
gmt plot basement_high.txt -W1p,0/0/150,- -V
echo 119.75 23.90 10,1,0/0/150 MC Basement | gmt text -F+f+j -V
echo 119.75 23.83 10,1,0/0/150 MC High | gmt text -F+f+j -V
REM ***profile - track & project***
echo 120.35 24.12 12,1 MC A | gmt text -F+f+j -N -V
echo 121.45 23.75 12,1 MC A'| gmt text -F+f+j -N -V
gmt plot track.A -W2p,0 -V
REM ****plot symbol, title ****
echo 119.75 25.70 8,1 MC Depth | gmt text -F+f+j -N -V
echo 120.00 25.70 8,1 MC 0-15 | gmt text -F+f+j -N -V
echo 120.25 25.70 8,1 MC 15-30 | gmt text -F+f+j -N -V
echo 120.50 25.70 8,1 MC 30-50 | gmt text -F+f+j -N -V
echo 120.74 25.70 8,1 MC \076 50 | gmt text -F+f+j -N -V
echo 120.90 25.70 8,1 MC km | gmt text -F+f+j -N -V
echo 121.15 25.80 10,1 LM Projection width\072 | gmt text -F+f+j -N -V
echo 121.15 25.69 10,1 LM 20 km (each side) | gmt text -F+f+j -N -V
echo 120.75 25.95 16,1 MC Background seismicities of Taiwan | gmt text -F+f+j -N -V
echo 120.00 25.78 | gmt plot -Sc.08i -W0.5p,255/106/106 -N -V
echo 120.25 25.78 | gmt plot -Sc.08i -W0.5p,250/250/0 -N -V
echo 120.50 25.78 | gmt plot -Sc.08i -W0.5p,50/250/100 -N -V
echo 120.74 25.78 | gmt plot -Sc.08i -W0.5p,0/191/255 -N -V
gmt colorbar -Ctopo_gray.cpt -Dx4.5c/-1c/8c/0.3c -By+l"m" -Bxa1000f500+l"Elevation" -V
echo Active faults (CGS, 2010) and Lungchung fault
gmt plot ActiveFaults_CGS2010_WGS84.txt -W0.8p,darkred -V
:: gmt plot LCNF.txt -W1.5p,black,- -V
gmt meca 19962024eqnew.txt -Sa0.4c+f6p -C+s1.5p -L -Fa2p/cc -Fe255/0/0 -V
gawk "{print $1, $2, $7*0.022}" 19962024eqnew.txt | gmt plot -Sa+gred -G0/0/0 -V
REM *****Plot strenghth profile based on seismicity*****
REM Line AA' profile & eq projection
gmt plot trackA.xyp -R0/118/0/3100 -JX3.5i/1.6i -G127.5 -W1p -Bxa20f10+ukm -Bya1000f500+l"Elev (m)" -BWNbr+t"Profile AA'" -X6i -Y5i -V
echo 5 2700 14 LB W | gmt text -F+f+j -N -V
echo 110 2700 14 LB E | gmt text -F+f+j -N -V
gmt coupe CWAEQs_1995-2022.txt -JX3.5i/1.5i -R0/118/-50/0 -Bxa20f10+l"Distance (Km)" -Bya10f5+l"Depth (km)" -BWStr -W0.5p ^
-Fsc0.04c/0 -Aa120.40/24.10/121.45/23.65/90/20/0/50 -Y-1.5i -V
REM ******Plot focal mechanisms under the profile******
gmt plot trackA.xyp -R0/118/0/3100 -JX3.5i/1.6i -G0 -W1p -Bxa20f10+ukm -Bya1000f500+l"Elev (m)" -BWNbr+t"Profile AA'" -Y-3i -V
echo 5 2700 14 LB W | gmt text -F+f+j -N -V
echo 110 2700 14 LB E | gmt text -F+f+j -N -V
gawk "{ print $1, $2, $3, $4, $5, $6, $7}" 19962024eqnew.txt > 19962024eqnew.gmt
gmt coupe 19962024eqnew.gmt -JX3.5i/1.5i -R0/118/-50/0 -Sa0.6c -Aa120.40/24.10/121.45/23.65/90/20/0/50 -G237/28/36 -E255 ^
-Bxa20f10+l"Distance (Km)" -Bya10f5+l"Depth (km)" -BWStr -W0.5p -Y-1.5i -V
gmt end
del *.nc *.conf *.history
pause
```
## 目前的圖

${X_{\alpha}} = {{{\lambda}h} \over {2 {\beta}{cos{\theta}}}}$
${{\Delta}{\phi}} = {-4{\pi} \over {\lambda}}{({R_M}-{R_S})}+{\Delta}{\phi}_{error}$
${{{\Delta}{\phi}}={2\pi \over \lambda}{\delta}}$
${({\gamma}+{\delta})^2}={\gamma}^2 + B^2 - 2{\gamma}Bcos({90^{\circ}}-{\theta}+{\alpha})$
${({\gamma}+{\delta})^2}={\gamma}^2 + B^2 - 2{\gamma}Bsin({\theta}-{\alpha})$
$sin({\theta}-{\alpha})={{{({\gamma}+{\delta})^2}-{\gamma}^2 - B^2} \over {2{\gamma}B}}$
$sin({\theta}-{\alpha}) = {{\delta} \over B}+{{{\delta}^2-B^2} \over {2{\gamma}B}}$
$Z(y) = h - {\gamma}cos{\theta}$
${Z(y)}=h-{{{{{\lambda}{\phi} \over 4{\pi}})^2}-B^2} \over 2B{sin}{{{\theta}-{\alpha}}-{({\lambda}{\phi} /4{\pi})}}$
${Z(y)}=h- {{(\frac {{\lambda}{\phi}}{4\pi})^2 - B^2} \over2(Bsin({\theta}-{\alpha})-\frac {{\lambda}{\phi}}{4\pi})}cos{\theta}$
${SCR} = {s^2 \over c^2}$
w