CPP-Library/RadiationPattern.hpp at master · shuleyu/CPP-Library · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#ifndef ASU_RADIATIONPATTERN
#define ASU_RADIATIONPATTERN

#include<iostream>
#include<cmath>
#include<string>

/*************************************************************
 * This C++ template calculate the radiation pattern amplitude
 * for given fault orientation and ray path orientation on
 * the lower hamisphere.
 *
 * input(s):
 * const double &strike    ----  Falut strike. (in deg, 0~360)
 * const double &dip       ----  Falut dip. (in deg, 0~90)
 * const double &rake      ----  Falut rake.  (in deg, -180~180)
 * const double &az        ----  Ray path azimuth. (in deg)
 * const double &take_off  ----  Ray path take-off angle. (in deg, 0~90)
 * const string &Cmp   ----  Chosen component. ("P","SV" or "SH").
 *
 * return(s):
 * double ans  ----  The radiation pattern intensity.
 *
 * Shule Yu
 * Jan 23 2018
 *
 * Key words: radiation pattern.
 *
 * Reference:
 *     Aki & Richard 2002, Page 108~109.
*************************************************************/

double RadiationPattern(const double &strike, const double &dip, const double &rake,
                        const double &az, const double &take_off,const std::string &Cmp){

    // Convert parameters into rad.
    double s,d,r,a,t;
    s=strike*M_PI/180;
    d=dip*M_PI/180;
    r=rake*M_PI/180;
    a=az*M_PI/180;
    t=take_off*M_PI/180;

    // Calculate.
    double ans;
    if (Cmp=="P"){

        ans=cos(r)*sin(d)*pow(sin(t),2)*sin(2*(a-s))-
            cos(r)*cos(d)*sin(2*t)*cos(a-s)+
            sin(r)*sin(2*d)*(pow(cos(t),2)-pow(sin(t),2)*pow(sin(a-s),2))+
            sin(r)*cos(2*d)*sin(2*t)*sin(a-s);

        return ans;
    }

    if (Cmp=="SV"){

        ans=sin(r)*cos(2*d)*cos(2*t)*sin(a-s)-
            cos(r)*cos(d)*cos(2*t)*cos(a-s)+
            0.5*cos(r)*sin(d)*sin(2*t)*sin(2*(a-s))-
            0.5*sin(r)*sin(2*d)*sin(2*t)*(1+pow(sin(a-s),2));

        return ans;
    }

    if (Cmp=="SH"){

        ans=cos(r)*cos(d)*cos(t)*sin(a-s)+
            cos(r)*sin(d)*sin(t)*cos(2*(a-s))+
            sin(r)*cos(2*d)*cos(t)*cos(a-s)-
            0.5*sin(r)*sin(2*d)*sin(t)*sin(2*(a-s));

        return ans;
    }

    std::cerr <<  "Error in " << __func__ << ": Component error (must be \"P\",\"SV\" or \"SH\") ... " << std::endl;

    return -1;
}

#endif