NINA/two_comp_field_decay_mod.cpp at master · ignotur/NINA · 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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#include <cmath>
#include "stars.h"

using namespace std;

//double int_xi (double, double, double, double, double, double);

//-----------------------------------------------------------------
// Период вращения нейтронной звезды. Формулы для периода с учётом
// убывания магнитного поля получены в основной работе методом
// интегрирования дифференциального уравнения для электромагнитного
// дипольного излучения и токовых потерь

double NeutronStar::get_P(double t/*, parametrs_B * param*/) {
    double P_res, I, tmp, tau_Hall, tau_ohm, t_step, I_dp, k_1, k_2, k_3, k_4, kx_1, kx_2, kx_3, kx_4, xi_res, diff_t;
    double param, intpart, fracpart, perm;
    t = t - tau;

    P_res  = P*P;


    if (t <= 3.5e5)												{
        param = t/1.e4;
        fracpart = modf(param, &intpart);
        t_step = 1.e4;

        for (int i=0; i < intpart; i++) 							{
            k_1 = t_step*3.2e7/2./1.e39*pow(get_B(tau+i*t_step), 2.);
            k_2 = t_step*3.2e7/2./1.e39*pow(get_B(tau+(i+0.5)*t_step), 2.);
            k_3 = t_step*3.2e7/2./1.e39*pow(get_B(tau+(i+0.5)*t_step), 2.);
            k_4 = t_step*3.2e7/2./1.e39*pow(get_B(tau+(i+1.0)*t_step), 2.);

            P_res += 0.1666667 * (k_1 + 2*k_2 + 2*k_3 + k_4);
        }

        t_step = fracpart*t_step;
        perm = intpart*t_step;

        k_1 = t_step*3.2e7/2./1.e39*pow(get_B(perm+tau), 2.);
        k_2 = t_step*3.2e7/2./1.e39*pow(get_B(perm+tau+0.5*t_step), 2.);
        k_3 = t_step*3.2e7/2./1.e39*pow(get_B(perm+tau+0.5*t_step), 2.);
        k_4 = t_step*3.2e7/2./1.e39*pow(get_B(perm+tau+1.0*t_step), 2.);

        P_res += 0.1666667 * (k_1 + 2*k_2 + 2*k_3 + k_4);
    } else if (t > 3.5e5)											{
        t_step = 1.e4;

        for (int i=0; i < 35; i++) 						{
            k_1 = t_step*3.2e7/2./1.e39*pow(get_B(tau+i*t_step), 2.);
            k_2 = t_step*3.2e7/2./1.e39*pow(get_B(tau+(i+0.5)*t_step), 2.);
            k_3 = t_step*3.2e7/2./1.e39*pow(get_B(tau+(i+0.5)*t_step), 2.);
            k_4 = t_step*3.2e7/2./1.e39*pow(get_B(tau+(i+1.0)*t_step), 2.);
            P_res += 0.1666667 * (k_1 + 2*k_2 + 2*k_3 + k_4);
        }

        diff_t = t - 3.5e5;
        P_res += 2./1.e39 * pow(get_B(3.7e5 + tau),2.)*diff_t*3.2e7;
    }


    return sqrt(P_res);
}

//-------------------------------------------------------------------


double NeutronStar::get_incl(double t) {
    double res, I_dp;

    return I_dp;
}

//-------------------------------------------------------------------
// Первая производная периода. Формулы взяты из исходного дифферен-
// циального уравнения

double NeutronStar::get_dot_P (double t/*, parametrs_B * param*/) {
    double res, I, tau_Hall, tau_ohm, I_dp;

    res = pow(get_B(t), 2.)/2e39/get_P(t);

    return res;
}
//-------------------------------------------------------------------


//------------------------------------------------------------------
// Падение магнитного поля в соотвествии со статьёй Pons, 2009
// Формула получена из личной переписки

double NeutronStar::get_B (double T/*, parametrs_B * param*/) {
    double res_B, B_min;
    double tau_Hall   = paramet_B->get_tau_Hall();
    double tau_ohm = paramet_B->get_tau_ohm();
    T = T - tau;

    //	if (T > 350e3)
    //		res_B=B*pow(pow(0.034*350e3/1e4, 1.17) + 0.84, -1.);
    //	else if (T <= 350e3 && T > 50e3)
    //		res_B=B*pow(pow(0.034*(T-50e3)/1e4, 1.17) + 0.84, -1.);
    //	else if (T < 50e3)
    //		res_B = B;


    res_B = B * exp(-T/5.e6);

    return res_B;
}
//------------------------------------------------------------------

parametrs_B::parametrs_B (ifstream * in) {
    *in>>tau_Hall;
    *in>>tau_ohm;
}
/*
parametrs_B::parametrs_B (parametrs_B * param)	{
alpha   = param->get_alpha();
tau_ohm = param->get_tau_ohm();
}
*/
void parametrs_B::print_description (ostream * out) {
    *out<<"//----------------------------------------------------------//"<<endl;
    *out<<"// Используется модель убывания магнитного поля с двумя эк- //"<<endl;
    *out<<"// споненциальными компонентами,   tau_Hall и вторая tau_ohm//"<<endl;
    *out<<"// так же учитывается убывание угла между магнитным полюсом //"<<endl;
    *out<<"// и полюсом вращения.                                      //"<<endl;
    *out<<"//----------------------------------------------------------//"<<endl;
}

void parametrs_B::print_parametrs   (ostream * out) {
    *out<<"//           Параметры модели убывания поля                 //"<<endl;
    *out<<"//----------------------------------------------------------//"<<endl;
    *out<<"// tau_Hall   - "<<tau_Hall<<endl;
    *out<<"// tau_ohm - "<<tau_ohm<<endl;
    *out<<"//----------------------------------------------------------//"<<endl;
}

void parametrs_B::print_short       (ostream * out) {
    *out<<"tau_Hall   - "<<tau_Hall<<endl;
    *out<<"tau_ohm - "<<tau_ohm<<endl;
}

double parametrs_B::get_tau_Hall (void) {
    return tau_Hall;
}

double parametrs_B::get_tau_ohm (void) {
    return tau_ohm;
}