辐射花样的计算与震源球的绘制

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14159265358979323846
#define DEG2RAD PI/180.0

int main (int argc, char *argv[]) {
    int i, j;
    float m[3][3];

    if (argc != 7 && argc != 4) {
        fprintf(stderr,"Usage: %s mrr mtt mff mrt mrf mtf\n", argv[0]);
        fprintf(stderr,"   Or: %s strike dip rake\n", argv[0]);
        exit(1);
    }

    if (argc == 7) {    // moment tensor
        sscanf(argv[1], "%f", &m[2][2]);    // mrr -> mzz
        sscanf(argv[2], "%f", &m[0][0]);    // mtt -> mxx
        sscanf(argv[3], "%f", &m[1][1]);    // mff -> myy

        sscanf(argv[4], "%f", &m[2][0]);    // mrt -> mzx
        m[0][2] = m[2][0];

        sscanf(argv[5], "%f", &m[2][1]);    // mrf -> -Mzy
        m[2][1] = -m[2][1];
        m[1][2] = m[2][1];

        sscanf(argv[6], "%f", &m[0][1]);    // mtf -> -Mxy
        m[0][1] = -m[0][1];
        m[1][0] = m[0][1];
    } else if (argc == 4) { // strike, dip, rake
        float strike, dip, rake;
        sscanf(argv[1], "%f", &strike);
        sscanf(argv[2], "%f", &dip   );
        sscanf(argv[3], "%f", &rake  );
        strike *= DEG2RAD;
        rake   *= DEG2RAD;
        dip    *= DEG2RAD;
        m[0][0] = - sin(dip)*cos(rake)*sin(2*strike)
                - sin(2*dip)*sin(rake)*sin(strike)*sin(strike);
        m[0][1] = sin(dip)*cos(rake)*cos(2*strike)
                + 0.5*sin(2*dip)*sin(rake)*sin(2*strike);
        m[0][2] = -cos(dip)*cos(rake)*cos(strike)
                - cos(2*dip)*sin(rake)*sin(strike);
        m[1][1] = sin(dip)*cos(rake)*sin(2*strike)
                - sin(2*dip)*sin(rake)*cos(strike)*cos(strike);
        m[1][2] = -cos(dip)*cos(rake)*sin(strike)
                + cos(2*dip)*sin(rake)*cos(strike);
        m[2][2] = sin(2*dip)*sin(rake);
        m[1][0] = m[0][1];
        m[2][0] = m[0][2];
        m[2][1] = m[1][2];
    }

    fprintf(stdout,"    / %6.3f %6.3f %6.3f \\ \n", m[0][0], m[0][1], m[0][2]);
    fprintf(stdout,"M = | %6.3f %6.3f %6.3f | \n",  m[1][0], m[1][1], m[1][2]);
    fprintf(stdout,"    \\ %6.3f %6.3f %6.3f / \n", m[2][0], m[2][1], m[2][2]);

    FILE *fop;
    fop = fopen("pattern.dat", "wb");
    double az, theta;
    float p[3];         // 离源矢量
    for (i=0; i<3600; i++)
      for (j=0; j<=900; j++) {
        az = (double)i / 10.0 * DEG2RAD;     // 方位角
        theta = (double)j/10.0 * DEG2RAD;    // 离源角，仅计算下半球
        p[0] = (float)(sin(theta)*cos(az));
        p[1] = (float)(sin(theta)*sin(az));
        p[2] = (float)(cos(theta));

        int k, l;
        float amp = 0.0;
        for (k=0; k<=2; k++)
          for (l=0; l<=2; l++){
            amp += p[k]*m[k][l]*p[l];
        }
        fwrite(&amp, sizeof(float), 1, fop);
    }
    fclose(fop);

    return 0;
}

#!/bin/bash
R=0/360/-90/0
J=S0/-90/15c
B=a30f10
name=pattern
PS=${name}.ps

gmt set MAP_FRAME_TYPE=plain
gmt set FORMAT_GEO_MAP=+D
gmt xyz2grd ${name}.dat -G${name}.nc -I6m/6m -R$R -ZLBxf
gmt grd2cpt ${name}.nc -Cpolar -E100 > ${name}.cpt
gmt psxy -R$R -J$J -T -K -P > $PS
gmt grdimage ${name}.nc -R$R -J$J -C${name}.cpt -B$B -BN -K -O >> $PS
gmt grdcontour ${name}.nc -R$R -J$J -L-0.001/0.001 -C1 -K -O -W2p >> $PS
gmt psxy -R$R -J$J -T -O >> $PS
rm gmt.conf gmt.history ${name}.cpt ${name}.nc 

#!/bin/bash
gmt psmeca -R0/250/-90/90 -JQ22c -Sm7c -W1p -Bx60 -By30 <<EOF> gmt_meca.ps
121.71 -12.01 29 -1.35 5.41 -4.06 -3.21 -3.58 -0.74 24 X Y 201304191751A
EOF

#!/bin/bash
gmt psmeca -R0/250/-90/90 -JQ22c -Sd7c -W1p -Bx60 -By30 <<EOF> gmt_meca.ps
121.71 -12.01 29 -1.35 5.41 -4.06 -3.21 -3.58 -0.74 24 X Y 201304191751A
EOF