NINA/Geppert_field_decay.cpp at master · Ingwar/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
#include <cmath>
#include "stars.h"

using namespace std;

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

double NeutronStar::get_P(double t/*, parametrs_B * param*/) {
    double P_res, I, tmp, alpha, tau_ohm;
    I = 2./5. * M*M_sol*pow(R,2);
    t = t - tau;

    alpha   = paramet_B->get_alpha();
    tau_ohm = paramet_B->get_tau_ohm();

    double tau_mu = 3.*pow(light_velocity,3)*I/2./pow(B,2)/pow(R,6);

    if (t<=1e6)	{
        tmp = alpha*pow(e, -t/tau_ohm);
        P_res = sqrt(P*P - 8*3.2e7*pi*pi*tau_ohm/tau_mu*(log(abs(tmp-alpha-1))/alpha/alpha-(alpha+1)/(tmp*alpha*alpha-alpha*alpha*alpha-alpha*alpha)) + 8*3.2e7*pi*pi*tau_ohm/tau_mu/alpha/alpha*(alpha+1));
    } else {
        tmp = get_P(tau+1e6);
        P_res = sqrt(tmp*tmp + 16*pi*pi*pow(R, 6)*pow(get_B(t),2)/3./I/pow(light_velocity, 3)*t/lsec);
    }

    return P_res;
}

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

double NeutronStar::get_incl (double t) {
    return i_incl;
}

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

double NeutronStar::get_dot_P (double t/*, parametrs_B * param*/) {
    double res, I, alpha, tau_ohm;
    alpha   = paramet_B->get_alpha();
    tau_ohm = paramet_B->get_tau_ohm();
    t = t - tau;
    I = 2./5. * M*M_sol*pow(R,2);
    double tau_mu = 3.*pow(light_velocity,3)*I/2./pow(B,2)/pow(R,6);
    /*	if (t>1e6)
    	res = pow(e, -2.*1e6/tau_ohm)/get_P(t+tau)/tau_mu/pow(1.+alpha*(1-pow(e, -1e6/tau_ohm)),2);
    	else
    	res = pow(e, -2.*t/tau_ohm)/get_P(t+tau)/tau_mu/pow(1.+alpha*(1-pow(e, -t/tau_ohm)),2);
    */

    if (t>1e6) {
        res = 4*pi*pi/tau_mu*pow(e, -2.*1e6/tau_ohm)/get_P(t+tau)/pow(1.+alpha*(1.-pow(e,-1e6/tau_ohm)), 2);
    } else {
        res = 4*pi*pi/tau_mu*pow(e, -2*t/tau_ohm)/get_P(t+tau)/pow(1.+alpha*(1.-pow(e,-t/tau_ohm)), 2);
    }


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


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

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

    if (T<1e6) {
        res_B = B * pow(e, -T/tau_ohm) / (1.+alpha*(1.-pow(e, -T/tau_ohm)));
    } else {
        res_B = B * pow(e, -1e6/tau_ohm) / (1.+alpha*(1.-pow(e, -1e6/tau_ohm)));
    }

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

parametrs_B::parametrs_B (ifstream * in) {
    *in>>alpha;
    *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<<"// работе Pons, et al. 2007. Первый миллион лет идёт экспо- //"<<endl;
    *out<<"// ненциальное убывание напряжённости магнитного поля, затем//"<<endl;
    *out<<"// напряжённость поля не меняется в ходе эволюции.          //"<<endl;
    *out<<"//----------------------------------------------------------//"<<endl;
}

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

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

double parametrs_B::get_alpha (void) {
    return alpha;
}

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