151 lines
6.2 KiB
C++
Executable File
151 lines
6.2 KiB
C++
Executable File
/*Calcule par interpolation linéaire des valeurs de flux stockées dans les fichiers de données "data_ord_p.res" les valeurs obtenues par le modèle pour les mêmes fréquences que celles mesurées dans "photandseds.tbl", le SED du quasar étudié, de sorte à pouvoir comparer les données et choisir le meilleur fit.*/
|
|
|
|
#include "functions.h"
|
|
#include<math.h>
|
|
#include<fstream>
|
|
#include<iostream>
|
|
#include<stdlib.h>
|
|
#include<time.h>
|
|
using namespace std;
|
|
|
|
const double e = 1.602e-19, me = 9.109e-31, eps0 = 1e-9/(36.*M_PI), c = 299792458;
|
|
double B1 = 1e-4, gam = 1e4, y_p1=pow(10,4.5), phi1 = M_PI/4., pitch = M_PI/4., kappa=1., p1=2.0, y_bulk1=4.0;
|
|
double gam2 = gam*gam;
|
|
double zeta = e*me/(24.*pow(M_PI,4)*eps0*c);
|
|
|
|
double z = 0.77, D_l = 3407.27e6*3.085e16; //paramètres redshift et distance lumière pour 3C175
|
|
|
|
/*
|
|
* Cas de l'électron isolé
|
|
*/
|
|
double beta(double y){ //calcul du facteur beta (v/c) en fonction du facteur de lorentz y
|
|
return sqrt(y*y-1.)/y;
|
|
}
|
|
|
|
double a(double B,double alpha, double y){ //calcul du rayon de courbure de la trajectoire de l'électron de pitch angle alpha et de facteur de lorentz y dans un champ magnétique d'intensité B
|
|
return sqrt(y*y-1)*me*c/(e*B*sin(alpha));
|
|
}
|
|
|
|
double thetagamma(double theta, double y){ //facteur provenant du développement limité du Flux en theta
|
|
return sqrt(1+y*y*theta*theta);
|
|
}
|
|
|
|
double eta(double w, double B, double alpha, double thetgam, double y){ //facteur provenant d'un changement de variable pour le calcul du flux
|
|
return (w*a(B,alpha,y)*pow(thetgam,3))/(3.*c*pow(y,3));
|
|
}
|
|
|
|
double L_isol(double w, double B, double phi, double alpha, double y){ //Cacul de l'émissivité pour un électron isolé
|
|
double theta = phi-alpha;
|
|
double thetgam = thetagamma(theta,y);
|
|
double eta_w = eta(w,B,alpha,thetgam,y);
|
|
double bet = beta(y);
|
|
|
|
double factor = (zeta*bet*bet*w*w*pow(thetgam,4))/(B*pow(y,3)*sin(alpha)*sin(alpha));
|
|
|
|
double K23 = Bessel_K23(eta_w), K13 = Bessel_K13(eta_w);
|
|
return factor*(K23*K23+(thetgam*thetgam-1)/(thetgam*thetgam)*K13*K13);
|
|
}
|
|
|
|
|
|
/*
|
|
* Pour une population d'électrons
|
|
*/
|
|
double pas_adapte(int d, double x, double x1, double x2, double y1, double y2){ //permet d'adapter le nombre de points d'intégration (donc le pas) selon une loi puissance négative
|
|
double a = (y1-y2)/(1./pow(x1,d)-1./pow(x2,d));
|
|
return a*(1./pow(x,d)-1./pow(x2,d))+y2;
|
|
}
|
|
|
|
|
|
double fonction(double w, double B, double phi, double alpha, double y, double p){ //expression de la fonction à intégrer selon alpha et gamma pour obtenir l'émissivité
|
|
return 0.5*sin(alpha)*kappa*pow(y,-p)*L_isol(w,B,phi,alpha,y);
|
|
}
|
|
|
|
double fy(double* tab, double y){ //fonction(w,phi,alpha,y) mais dont seul 'y' est variable, pour l'intégration selon 'y'
|
|
return fonction(tab[0],tab[1],tab[2],tab[4],y,tab[5]);
|
|
}
|
|
|
|
double fa(double* tab, double alpha){ //fonction(w,phi,alpha,y) mais dont seul 'alpha' est variable, pour l'intégration selon 'alpha'
|
|
double w = tab[0], y_p = tab[3];
|
|
tab[4] = alpha;
|
|
int n = (int)pas_adapte(6,log10(w),7.,23.,1e6,50.);
|
|
return GaussIntegration2(1.,y_p,n,fy,tab);
|
|
}
|
|
|
|
double L_pop(double w, double B, double phi, double y_p, double p, double* tab){ //calcul de l'émissivité pour une population d'électrons
|
|
double dalpha = 1e-3/sqrt((w/(2*M_PI))*1e-10);
|
|
double alpha1 = phi-dalpha, alpha2 = phi+dalpha;
|
|
int n = 50;//(int)pas_adapte(5,log10(w),7.,23.,1000.,50.);
|
|
tab[0] = w ; tab[1] = B ; tab[2] = phi ; tab[3] = y_p; tab[4] = alpha1; tab[5] = p;
|
|
return GaussIntegration2(alpha1,alpha2,n,fa,tab);
|
|
}
|
|
|
|
/*
|
|
* Boost relativiste
|
|
*/
|
|
|
|
double DopplerFactor(double y, double phi){
|
|
double bet = beta(y);
|
|
return 1./(y*(1-bet*cos(phi)));
|
|
}
|
|
|
|
double phi_rest(double phi, double y){ //calcul de phi dans le référentiel dans le référentiel au repos d'un électron de facteur de lorentz y dans le référentiel de l'observateur
|
|
double d = DopplerFactor(y,phi);
|
|
return acos((cos(phi)-beta(y))/(1-beta(y)*cos(phi)));
|
|
//return asin(d*sin(phi));
|
|
}
|
|
|
|
int main(){
|
|
time_t start = time(nullptr);
|
|
fstream data("Flux_modele.res",ios::out);
|
|
double wi = 1e7, wf = 1e15;
|
|
double phi_i = 1e-3, phi_f = M_PI/2.;
|
|
int i,j, n=100;
|
|
double h = (log10(wf)-log10(wi))/n, g = (phi_f-phi_i)/n;
|
|
double* m = (double*)malloc(6*sizeof(double)); //pointeurs dont on va avoir besoin dans L_pop, fy(*m,y) et fa(*m,alpha) pour se passer les arguments fixes
|
|
double* l = (double*)malloc(6*sizeof(double));
|
|
|
|
data << phi1;
|
|
|
|
double *p = (double*)malloc(5*sizeof(double)); //valeurs de p testées
|
|
for(j=0;j<5;j++){
|
|
p[j] = 2.75+0.125*j;
|
|
data << " " << p[j];
|
|
}
|
|
data << endl;
|
|
|
|
double d,w1=0,w2=0,f,u1,u2;
|
|
double mem[8] = {0.,0.,0.,0.,0.,0.,0.,0.}, fromf1[9] = {0.,0.,0.,0.,0.,0.,0.,0.,0.}, fromf2[9] = {1.,0.,0.,0.,0.,0.,0.,0.,0.}; //mem stocke un point interpolé dans l'intervalle [fromf1;fromf2]
|
|
fstream fich("photandseds.tbl",ios::in);
|
|
for(j=0;j<4;j++){
|
|
fich.ignore(1000,'\n');
|
|
}
|
|
while(fich >> w1 >> f >> u1 >> u2){ //on récupère les abscisses des points de "photandseds.tbl"
|
|
data << w1;
|
|
if(w1 != w2){
|
|
fstream forinterpol("data_ord_p_3.res",ios::in);
|
|
forinterpol.ignore(1000,'\n');
|
|
while(forinterpol >> fromf2[0] >> fromf2[1] >> fromf2[2] >> fromf2[3] >> fromf2[4] >> fromf2[5]){ //on récupère les valeurs fromf1 et fromf2 du set de données "data_ord_p.res"
|
|
if(fromf2[0] > w1)
|
|
break;
|
|
for(j=0;j<6;j++)
|
|
fromf1[j] = fromf2[j];
|
|
}
|
|
forinterpol.close();
|
|
for(j=0;j<5;j++){ //on effectue l'interpolation linéaire des points de "data_ord_p.res" aux abscisses des points de "photandseds.tbl"
|
|
//d = DopplerFactor(y_bulk1,phi1);
|
|
mem[j] = fromf1[j+1]+(fromf2[j+1]-fromf1[j+1])*(w1-fromf1[0])/(fromf2[0]-fromf1[0]);//pow(d,2.+(p[j]-1.)/2.)*L_pop(2.*M_PI*w1,B1,phi_rest(phi1,y_bulk1),y_p1,p[j],m)*pow(1.+z,(p[j]+1.)/2.)/(4.*M_PI*D_l*D_l);
|
|
w2 = w1;
|
|
}
|
|
}
|
|
for(j=0;j<5;j++)
|
|
data << " " << mem[j];
|
|
data << endl;
|
|
}
|
|
fich.close();
|
|
data.close();
|
|
free(m); free(l); free(p);
|
|
cout << "Elapsed time : " << difftime(time(nullptr),start) << "s." << endl;
|
|
system("python plotfromsed.py"); //tracé du fit avec python
|
|
return 0;
|
|
}
|