/*
Biblietoile - bibliothèque de routine pour calcul de structure stellaire
Copyright (C) 2002 Observatoire de Paris-Meudon
Ce programme est un logiciel libre ; vous pouvez le redistribuer et/ou le modifier conformément aux dispositions de la Licence Publique Générale GNU, telle que publiée par la Free Software Foundation ; version 2 de la licence, ou encore (à votre choix) toute version ultérieure.
Ce programme est distribué dans l'espoir qu'il sera utile, mais SANS AUCUNE GARANTIE ; sans même la garantie implicite de COMMERCIALISATION ou D'ADAPTATION A UN OBJET PARTICULIER. Pour plus de détail, voir la Licence Publique Générale GNU .
Vous devez avoir reçu un exemplaire de la Licence Publique Générale GNU en même temps que ce programme ; si ce n'est pas le cas, écrivez à la Free Software Foundation Inc., 675 Mass Ave, Cambridge, MA 02139, Etats-Unis.
*/
import java.util.*;
/** Bibliothèque de routines pour le calcul d'une structure stellaire. Utilisées notamment par les
classes Etoile.java et StellarApplet.java.
*/
public class Biblietoile {
/*+-------------------------------------+
| Déclaration des Constantes globales |
+-------------------------------------+*/
static final double C = 2.998e10; // vitesse de la lumiére en cm/s
static final double G = 6.67e-8; // constante de la gravitation en cm3/g/s2
static final double SIGMA = 5.67e-5; // constante de Bolzmann en erg/cm2/s/K4
static final double ASTRO_UNIT = 1.4e13; // unité astronomique en cm
static final double SOLEIL_MASSE = 1.99e33; // masse solaire en g
static final double SOLEIL_LUM = 3.83e33; // luminosité solaire en erg/sec
static final double SOLEIL_RAYON = 6.96e10; // rayon du soleil en cm
static final double SOLEIL_AGE = 4.6e9; // age du soleil en années
static final double SOLEIL_T_CENTRAL = 15.5e6; // Temp. centrale du soleil en K
static final double SOLEIL_T_SURFACE = 5800; // Temp. de surface du soleil en K
static final double SOLEIL_D_CENTRAL = 147.2; // densité centrale du soleil en g/cm3
static final double SOLEIL_X = 0.71; // X du soleil
static final double SOLEIL_Y = 0.27; // Y du soleil
static final double SOLEIL_Z = 0.02; // Z du soleil
static final double FRACTION_D_AZOTE = 0.3; // fraction d'azote dans les éléments lourds
static final double ENERGIE_H_EN_HE = 7e18; // production d'énergie en convert. 1g d'H en He
/*+------------------------------------+
| Déclaration des variables globales |
+------------------------------------+*/
double fitPrecision; // precision désirée au point de raccordement
double modelPrecision; // precision désirée pour le modéle (=10*precision)
int modIteration; // compteur d'itérations
int plotControl; // pour sortir les profils T, D, R, L : 0=pas de transfert des valeurs,1=oui
boolean plotMasseOuRayon; // Tracés en fonction de : false -> masse, true -> rayon
double pi; // 3.14159...
double twoTo14; // 2^.25
int modeleCalcule; // égal à 0 suf si les valeurs calculées sont les valeurs courantes
double x; // fraction de masse d'hydrogéne
double y; // fraction de masse d'hélium
double z; // fraction de masse des "métaux"
double zAzote; // fraction de masse d'azote
double density; // Densité (g/cm3)
double temperature; // Température (K)
double mass; // Masse de l'étoile
double radius; // Rayon à la surface
double luminosite; // Luminosité à la surface
double tSurface; // Température de surface
double tEffective; // température effective de l'étoile
double tCentral; // température centraler
double dCentral; // densité centrale
double tZero; // "vraie" température au centre
double dZero; // "vraie" densité au centre
double dPcorr; // correction de pression
double gradlTlP; // =dlnTdlnP mais en global
double gradlDlP; // =dlnDdlnP mais en global
double rayonCoeur; // pour récupérer la valuer calculée dans dehors pour extrapolation AR
double fitMassAuto; // masse au point de jonction calculé automatiquement
double minDeltaGlob; // le minimum maximorum pour recherche auto. du point de jonction
double fitMass; // masse au point de jonction
double fitMassFraction; // fraction de masse de l'étoile au point de jonction
double densiteSoleil; // densité solaire moyenne
double constGrrad; // constante sun_lum/sun_mass/(16*pi*c*G) // approx. 2.e-5 cgs
double constdPdM; // constante G/4/pi*(Sun_mass/Sun_radius^2)^2
double kOpacite; // Opacité (cm^2/g)
double kElectron; // Opacité due à la diffusion par les électrons
double kKramer; // Opacité de Kramers
double kHMoins; // Opacité de l'ion H-
double kMoleculaire; // Opacité moléculaire
double kRadiative; // Opacité radiative totale
double kConductive; // Opacité due à la conductivité électronique
double puissancepp; // puissance engendrée par p-p
double puissanceCNO; // puissance engendrée par CNO
double puissance3alpha; // puissance engendrée par 3 alpha
double puissance; // puissance engendrée par p-p + CNO
double dxdt; // taux de déplétion de l'hydrogéne
double dydt; // taux de déplétion de l'hélium
double pression; // pression totale, rayonnement+ions+électrons (dyne/cm^2)
double dLogPdLogT; // d log(p)/d log(T) à densité constante
double dLogPdLogD; // d log(p)/d log(density) à température constante
double pressionRad; // pression de rayonnement
double pressionIon; // pression due aux ions
double pressionElect; // pression électronique
double pressionGaz; // pression du gaz, ions+électrons
double pressionDegR; // pression des électrons dégénérés relativistes
double pressionDegNR; // pression des électrons dégénérés non-relativistes
double pressionDeg; // pression des électrons dégénérés
double pressionNon; // pression des électrons non-dégénérés, loi des gaz parfaits
double dLnTdLnPrad; // gradient radiatif de température
double dLnTdLnPad; // d logT/d logP à entropie constante ; gradient de T adiabatique
double maxParamStep[] = new double[5]; // Incrément/décrément max. pour chacun des paramétres
double corrections[] = new double[5]; // corrections calculées par correctionValeur
double maxCorrect[] = new double[5]; // corrections max. effectuées par correctionValeur
/*+------------------------------------+
| Déclaration des sorties graphiques |
+------------------------------------+*/
static final int keepSize = 50; // dimension max des tableaux graphiques
double M_array[] = new double[keepSize]; // profil de masse
double T_array[] = new double[keepSize]; // profil de température
double D_array[] = new double[keepSize]; // profil de densité
double R_array[] = new double[keepSize]; // profil du rayon
double L_array[] = new double[keepSize]; // profil de luminosité
double K_array[] = new double[keepSize]; // profil d'opacité
double E_array[] = new double[keepSize]; // profil de taux de production d'énergie
double P_array[] = new double[keepSize]; // profil de pression
Vector vecteurXin = new Vector(); // vecteur pour stocker les abscisses des profils produits par dedans
Vector vecteurTin = new Vector(); // profil de température
Vector vecteurDin = new Vector(); // profil de densité
Vector vecteurRin = new Vector(); // rayon ou profil de masse
Vector vecteurLin = new Vector(); // profil de luminosité
Vector vecteurXout = new Vector(); // comme ci-dessus mais pour dehors
Vector vecteurTout = new Vector();
Vector vecteurDout = new Vector();
Vector vecteurRout = new Vector();
Vector vecteurLout = new Vector();
/** Méthode pour convertir les nombres réels en chaîne de caractères avec choix du format d'affichage.
On entre le nombre, le nombre de décimales souhaitées, et un booléen à true pour avoir un exposant.
*/
public String arrondi(double nombre, int decimales, boolean scienti) {
double exposant = 0, mantisse = 0;
String nombreArrondi = null;
boolean signe = false;
if ( nombre != 0 ) {
if ( nombre < 0 ) {
nombre = -nombre;
signe = true;
}
if ( !scienti ) {
nombre = nombre * Math.pow(10.,(double)decimales);
nombre = (double)Math.round(nombre);
nombre = nombre / Math.pow(10.,(double)decimales);
nombreArrondi = Double.toString(nombre);
}
if ( scienti | (Math.log(nombre)/Math.log(10.)) >= 6. ) {
exposant = Math.floor(Math.log(nombre)/Math.log(10.));
mantisse = nombre/Math.pow(10.,exposant);
nombreArrondi = arrondi(mantisse,decimales,false)+"E"+Integer.toString((int)exposant);
}
if ( signe ) {
nombreArrondi = "-"+nombreArrondi;
}
}
else {
nombreArrondi = "0.";
}
return nombreArrondi;
}
/** Pour initialiser le calcul : composition chimique, quelques constantes utiles.
*/
public void initialiseModele(){
// Variable(s) locale(s)
// ---------------------
int i = 0;
double intermediaire = 0.; // variable de type double pour tests
pi = 4.0*Math.atan(1);
twoTo14 = Math.pow(2.,0.25); // 2^.25
x = SOLEIL_X; // fraction de masse d'hydrogène
y = SOLEIL_Y; // fraction de masse d'hélium
z = SOLEIL_Z; // fraction de masse des "métaux"
zAzote = z*FRACTION_D_AZOTE; // fraction de masse d'azote
density = SOLEIL_D_CENTRAL; // densité
temperature = SOLEIL_T_CENTRAL; // température
radius = Math.sqrt(luminosite)/Math.pow((tEffective/SOLEIL_T_SURFACE),2);
constGrrad = (SOLEIL_LUM/SOLEIL_MASSE)/(16.0*pi*C*G); // utilisé dans EvaluateEq
constdPdM = G/4.0/pi*Math.pow((SOLEIL_MASSE/Math.pow(SOLEIL_RAYON,2)),2); // utilisé dans EvaluateEq
densiteSoleil = SOLEIL_MASSE/(4.0*pi*SOLEIL_RAYON*SOLEIL_RAYON*SOLEIL_RAYON/3.0); // densité moyenne du Soleil
maxParamStep[0] = 0.20; // Pas maximal pour la température - au lieu de .05
maxParamStep[1] = 0.20; // Pas maximal pour la densité - au lieu de .15
maxParamStep[2] = 0.20; // Pas maximal pour le rayon - au lieu de .05
maxParamStep[3] = 0.20; // Pas maximal pour la luminosité
maxParamStep[4] = 0.20; // Pas maximal pour la masse
fitPrecision = 0.001; // précision souhaitée au point de jonction
modelPrecision = 10.*fitPrecision; // précision souhaitée pour le modèle
// corrections max. pour chaque paramètre dans correctionValeur
maxCorrect[0] = 0.03; // Teff
maxCorrect[1] = 0.10; // Lum
maxCorrect[2] = 0.02; // Tcen - au lieu de .01
maxCorrect[3] = 0.10; // Dcen
maxCorrect[4] = 0.00; // Masse
// nettoyage des tableaux
double M_array = 0.;
double T_array = 0.;
double D_array = 0.;
double R_array = 0.;
double L_array = 0.;
double K_array = 0.;
double E_array = 0.;
double P_array = 0.;
}
/*+--------------------------------------------+
| Méthodes de calcul des processus physiques |
+--------------------------------------------+*/
/** Calcule la production d'énergie par la chaîne p-p et le cycle CNO.
Expression de KW. Accord avec BP sauf pour les facteurs d'écrantage
d'ordre élevé.
Entrée : x,y,z,zAzote,densite,temperature
Sortie : puissancepp, puissanceCNO, puissance (erg/g/sec),
dxdt (g/g/sec) (taux de disparition de H)
*/
public void energieNucleaire(double xx,double yy,double zz,double zn,double d,double t){
// Variable(s) locale(s)
// ---------------------
double t6; // temp./10^6
double t6to13; // t6 à la puissance 1/3
double t6to23; // t6 à la puissance 2/3
double screenpp; // facteur d'écrantage pour pp
double ecranCNO; // facteur d'écrantage pour cno
double tempo;
puissancepp = 1.0e-30; // pas de production d'énergie à basse T
puissanceCNO = 1.0e-30;
puissance = 1.0e-30;
if (t > 1.0e6) {
t6 = t/1.0e6; // T6
t6to13 = Math.pow(t6,0.3333); // T6 à la puissance 1/3
t6to23 = t6to13*t6to13; // T6 à la puissance 2/3
tempo= (64.33-152.28/t6to13); // au lieu de 64.24, 152.313
if (tempo > -30.0) {
ecranCNO = 1+0.0027*t6to13-0.00778*t6to23-0.000149*t6;
if (ecranCNO < 1.0) {
ecranCNO = 1.0;
}
puissanceCNO = xx*zn*ecranCNO*d*Math.exp(tempo)/t6to23;
}
screenpp = 1+0.0123*t6to13+0.0109*t6to23+0.0009*t6;
puissancepp = xx*xx*screenpp*d*Math.exp(14.68-33.80/t6to13)/t6to23; // au lieu de 14.54,33.81
puissance = puissancepp+puissanceCNO;
}
dxdt = -1.667e-19*puissance; // taux de disparition de l'hydrogène
}
/** Production d'énergie par la réaction 3 alpha.
Expression de KW. Accord avec BP sauf pour les facteurs d'écrantage
d'ordre élevé.
Entrée : x,y,z,zAzote,densite,temperature
Sortie : puissance3alpha (erg/g/sec),
dydt (g/g/sec) (taux de disparition de He)
*/
public void energieNucleaireAlpha(double xx,double yy,double zz,double zn,double d,double t){
// Variable(s) locale(s)
// ---------------------
double t8;
puissance3alpha = 1.0e-30;
if(t > 1.0e8) {
t8 = t/1.0e8;
puissance3alpha = 5.09e11*(d*d)*Math.pow((yy/t8),3)*Math.exp(-44.027/t8);
}
dydt = -1.7e-18*puissance3alpha; // taux de disparition de l'hélium
}
/** Calcul des opacités.
Entrée : x,y,z,densite,temperature
Sortie : opacité (cm^2/g) et ses différentes composantes
*/
public void opacite(double xx,double yy,double zz,double d,double t){
// Variable(s) locale(s)
// ---------------------
double zAverage;
// Constantes(s) locale(s)
// -----------------------
final double AA = 2; // d'après BP au lieu de 6.0 selon Kramers
final double BB = 1; // d'après BP au lieu de .001 H-
final double CC = 1; // d'après BP au lieu de .001 molecular
zAverage = xx+2*yy+8*zz; // au lieu de 4*yy ?
/* Calcul de l'opacité due à la diffusion électronique à haute T avec correction
de densité - moyennement précise (terme dominant selon KW, correction selon
BP d'après AC)
*/
kElectron = 0.2*(1+xx)/((1+2.7e11*d/t/t)*(1+Math.pow((t/4.5e8),0.86)));
/* Opacité de Kramer : décrit approximativement l'opacité des transitions
libre-libre, lié-libre, lié-lié.
*/
// kKramer = (13.7/AA)*(1+xx)*(0.0010+zz)*d*Math.pow((1.0e7/t),3.5) // formule de BP
kKramer = (4.56/AA)*(1+xx)*(0.0025+zz)*d*Math.pow((1.0e7/t),3.5); // formule de Stein
/* Opacités approximatives de H- et moléculaire
*/
if(t > 4.0e4) {
kHMoins = 65.0*Math.pow((zz*d),2)*Math.pow((4.0e4/3000.0),7.7); // valeur constante
kMoleculaire = 1.0e-5;
}
else {
kHMoins = 65.0*Math.pow((zz*d),2)*Math.pow((t/3000.0),7.7); // BP multiple ceci par BB
kMoleculaire = 0.1*zz; // BP multiple ceci par CC
}
// Opacité totale radiative
kRadiative = kMoleculaire + 1./(1./kHMoins+1./(kElectron+kKramer));
// Opacité par conductivité électronique, approximative
if(d < 1.0e-5) {
kConductive = 1.0e10;
}
else {
kConductive = zAverage*(2.6e-7)*Math.pow((t/d),2)*(1+Math.pow((d/2.0e6),0.66667));
}
// Opacité totale
kOpacite = 1./(1./kRadiative+1./kConductive);
}
/** Calcul de l'équation d'état et des grandeurs thermodynamiques. Expressions selon BP.
Entrée : x,y,z,densité,température
Sortie : pression, pressionGaz, pressionRad (dynes/cm^2),
dLogPdLogT, dLogPdLogD, dLnTdLnPad
L'ionisation totale est supposée, les électrons peuvent être partiellement
dégénérés et partiellement relativistes
*/
public void equationDEtat(double xx,double yy,double zz,double d,double t){
// Constantes(s) locale(s)
// -----------------------
final double radconst = 2.523e-15; // rad const a/3 (erg/CC/K^4)
final double ionconst = 0.8263e8; // const de Boltzmann/masse du proton (erg/g/K)
final double const1 = 1.0036e13; // KW -pour gaz d'e- dégénéré non-relativiste
final double const2 = 1.2435e15; // KW -pour gaz d'e- dégénéré relativiste
// Variable(s) locale(s)
// ---------------------
double poidsMolMoyenIon; // Poids moléculaire moyen des ions; on prend z de l'oxygène
double poidsMolMoyenElec; // Poids moléculaire moyen par électron libre, ionisation totale
double relativiste; // A quel point les e- sont-ils relativistes ?
double dLogPEdLogT; // d log(P elect)/d log(T) à densité constante
double dLogPEdLogD; // d log(P elect)/d log(densité) à T constante
double qt; // T*Cv, où Cv est la chaleur spécifique à V constant (par gramme)
double qr; // qt*(d log(T)/d log(densité)) à entropie constante
poidsMolMoyenIon = 1.0/(xx+yy/4.0+zz/16.0); // Poids moléculaire moyen des ions
poidsMolMoyenElec = 2.0/(1.0+xx); // Poids moléculaire moyen par électron libre
pressionRad = radconst*(Math.pow(t,4)); // pression du rayonnement
pressionIon = d*t*ionconst/poidsMolMoyenIon; // pression des ions
pressionDegR = const2*Math.pow((d/poidsMolMoyenElec),(4.0/3.0));
pressionDegNR = const1*Math.pow((d/poidsMolMoyenElec),(5.0/3.0));
pressionDeg = pressionDegNR/Math.pow((1.0+Math.pow((pressionDegNR/pressionDegR),2)),2);
pressionNon = d*t*ionconst/poidsMolMoyenElec;
pressionElect = pressionNon*Math.pow((1.0+Math.pow((pressionDeg/pressionNon),2)),2);
pressionGaz = pressionIon+pressionElect; // pression du gaz
pression = pressionGaz+pressionRad; // pression totale
// pression doit être de l'ordre de 10^17 dynes/cm2 au centre du soleil
relativiste = (5.0/3.0)*Math.pow((pressionDeg/pressionDegNR),2) + (4.0/3.0)*Math.pow((pressionDeg/pressionDegR),2);
dLogPEdLogT = Math.pow((pressionNon/pressionElect),2);
dLogPEdLogD = dLogPEdLogT+(1.0-dLogPEdLogT)*relativiste;
dLogPdLogT = (4.0*pressionRad+pressionIon + dLogPEdLogT*pressionElect)/pression;
dLogPdLogD = (pressionIon+dLogPEdLogD*pressionElect)/pression;
qt = (12.0*pressionRad+1.5*pressionIon+pressionElect*dLogPEdLogT/(relativiste-1.0))/d;
qr = pression*dLogPdLogT/d;
dLnTdLnPad = 1.0/(dLogPdLogD*qt/qr+dLogPdLogT); // DEL ADIABATIQUE
dLnTdLnPrad = constGrrad*kOpacite*pression/pressionRad; // DEL RAD*M/L
}
/** Résout les équations différentielles de la structure interne.
Entrée : tableau de valeurs initiales
Sortie : tableau des dérivées, d Params(i)/d Params(1)
Les paramètres sont : température, densité, rayon, luminosité, masse
*/
public Vecteur5 evalueEquations(Vecteur5 params){
Vecteur5 deriv = new Vecteur5();
// Variable(s) locale(s)
// ---------------------
double t; // température
double d; // densité
double r; // rayon
double l; // Luminosité
double m; // masse
double dLnTdLnP; // gradient de temperature dans l'étoile dLnt/dLnp
double dLnDdLnP; // gradient de densité dLnDensity/dLnp
double dPdM; // gradient de pression dLnp/d(M/Msun)
t = params.p0; // température
d = params.p1; // densité
r = params.p2; // rayon
l = params.p3; // Luminosité
m = params.p4; // masse
this.opacite(x,y,z,d,t); // on calcule d'abord l'opacité
this.equationDEtat(x,y,z,d,t); // calcul de l'équation d'état
this.energieNucleaire(x,y,z,zAzote,d,t); // calcul de la production d'énergie
dLnTdLnPrad = dLnTdLnPrad*l/m; // DEL RAD
// on prend le plus petit des gradients de T entre adiabatique et radiatif
dLnTdLnP = Math.min(dLnTdLnPad,dLnTdLnPrad);
dLnDdLnP = (1.-dLnTdLnP*dLogPdLogT)/dLogPdLogD; // gradient de densité
dPdM = -constdPdM*m/Math.pow(r,4); // gradient de pression
deriv.p0 = dLnTdLnP*dPdM*t/pression; // dT/dM
deriv.p1 = dLnDdLnP*dPdM*d/pression; // dD/dM
deriv.p2 = densiteSoleil/3.0/d/(r*r); // dr/dM
deriv.p3 = SOLEIL_MASSE/SOLEIL_LUM*puissance; // dL/dM
deriv.p4 = 1.; // dM/dM
return deriv;
}
/** Intégration d'un point de grille au suivant par méthode de Runge-Kutta à l'ordre 4 des grandeurs
passées en paramètre (de type Vecteur5).
Entrée : valeurs initiales
Sortie : valeurs après intégration sur un pas de grille
*/
public Vecteur5 integre1pas(Vecteur5 parameters) {
// Variable(s) locale(s)
// ---------------------
Vecteur5 derivs = new Vecteur5(); // derivées des paramètres
Vecteur5 newParameters = new Vecteur5(); // valeurs aux points intermédiaires RK
Vecteur5 k1 = new Vecteur5();
Vecteur5 k2 = new Vecteur5();
Vecteur5 k3 = new Vecteur5();
Vecteur5 k4 = new Vecteur5();
double pasInteg = 0; // pas d'intégration
derivs = evalueEquations(parameters); // résout les équations différentielles avec les paramètres
// Calcul du pas d'integration. Le pas d'intégration est le plus petit de ceux font varier les parametres de moins de maxparamstep en valeurs relatives
double maxinter = 0;
pasInteg = Math.abs(derivs.p0/maxParamStep[0]/parameters.p0);
pasInteg = Math.max(pasInteg,maxinter);
maxinter = pasInteg;
pasInteg = Math.abs(derivs.p1/maxParamStep[1]/parameters.p1);
pasInteg = Math.max(pasInteg,maxinter);
maxinter = pasInteg;
pasInteg = Math.abs(derivs.p2/maxParamStep[2]/parameters.p2);
pasInteg = Math.max(pasInteg,maxinter);
maxinter = pasInteg;
pasInteg = Math.abs(derivs.p3/maxParamStep[3]/parameters.p3);
pasInteg = Math.max(pasInteg,maxinter);
maxinter = pasInteg;
pasInteg = Math.abs(derivs.p4/maxParamStep[4]/parameters.p4);
pasInteg = Math.max(pasInteg,maxinter);
pasInteg = 1./pasInteg;
if (parameters.p4 > fitMass) {
pasInteg = -pasInteg; // pas négatif si masse > masse au point de jonction
}
if ((parameters.p4-fitMass)*(parameters.p4+pasInteg-fitMass) < 0.0) {
pasInteg = fitMass-parameters.p4; // pour le point le plus proche du point de jonction
}
// Calcul de k1
k1.p0 = pasInteg*derivs.p0;
k1.p1 = pasInteg*derivs.p1;
k1.p2 = pasInteg*derivs.p2;
k1.p3 = pasInteg*derivs.p3;
k1.p4 = pasInteg*derivs.p4;
// Calcul de k2
newParameters.p0 = parameters.p0 + k1.p0/2.;
newParameters.p1 = parameters.p1 + k1.p1/2.;
newParameters.p2 = parameters.p2 + k1.p2/2.;
newParameters.p3 = parameters.p3 + k1.p3/2.;
newParameters.p4 = parameters.p4 + k1.p4/2.;
derivs = evalueEquations(newParameters);
k2.p0 = pasInteg*derivs.p0;
k2.p1 = pasInteg*derivs.p1;
k2.p2 = pasInteg*derivs.p2;
k2.p3 = pasInteg*derivs.p3;
k2.p4 = pasInteg*derivs.p4;
// Calcul de k3
newParameters.p0 = parameters.p0 + k2.p0/2.;
newParameters.p1 = parameters.p1 + k2.p1/2.;
newParameters.p2 = parameters.p2 + k2.p2/2.;
newParameters.p3 = parameters.p3 + k2.p3/2.;
newParameters.p4 = parameters.p4 + k2.p4/2.;
derivs = evalueEquations(newParameters);
k3.p0 = pasInteg*derivs.p0;
k3.p1 = pasInteg*derivs.p1;
k3.p2 = pasInteg*derivs.p2;
k3.p3 = pasInteg*derivs.p3;
k3.p4 = pasInteg*derivs.p4;
// Calcul de k4
newParameters.p0 = parameters.p0 + k3.p0;
newParameters.p1 = parameters.p1 + k3.p1;
newParameters.p2 = parameters.p2 + k3.p2;
newParameters.p3 = parameters.p3 + k3.p3;
newParameters.p4 = parameters.p4 + k3.p4;
derivs = evalueEquations(newParameters);
k4.p0 = pasInteg*derivs.p0;
k4.p1 = pasInteg*derivs.p1;
k4.p2 = pasInteg*derivs.p2;
k4.p3 = pasInteg*derivs.p3;
k4.p4 = pasInteg*derivs.p4;
// calcul final
parameters.p0 = parameters.p0 + (k2.p0+k3.p0)/3.0 + (k1.p0+k4.p0)/6.0;
parameters.p1 = parameters.p1 + (k2.p1+k3.p1)/3.0 + (k1.p1+k4.p1)/6.0;
parameters.p2 = parameters.p2 + (k2.p2+k3.p2)/3.0 + (k1.p2+k4.p2)/6.0;
parameters.p3 = parameters.p3 + (k2.p3+k3.p3)/3.0 + (k1.p3+k4.p3)/6.0;
parameters.p4 = parameters.p4 + (k2.p4+k3.p4)/3.0 + (k1.p4+k4.p4)/6.0;
return parameters;
}
/** Intégration de la structure stellaire depuis le centre de l'étoile.
Entrée : température et densité au centre
Sortie : les 5 valeurs au point de jonction (T, densité, rayon, luminosité, masse)
*/
public Vecteur5 dehors(double tCent,double dCent,int error){
Vecteur5 core = new Vecteur5();
// Constantes(s) locale(s)
// -----------------------
final int maxiteration = 1000; // nombre maximum d'itérations
// Variable(s) locale(s)
// ---------------------
int iteration; // compteur d'itérations
double deviation; // écart de la masse au dernier point à la masse au point de jonction
double masseCoeur; // masse du coeur (en masse solaire)
double luminositeCoeur; // luminosité du coeur (en luminosité solaire)
double xx, yy;
int index;
this.energieNucleaire(x,y,z,zAzote,dCent,tCent); // calcul de la production d'énergie
masseCoeur = fitMass/10000.0; // masse du coeur à partir de laquelle on intègre
luminositeCoeur = SOLEIL_MASSE/SOLEIL_LUM*masseCoeur*puissance;
rayonCoeur = Math.pow((masseCoeur/dCent*densiteSoleil),(1.0/3.0));
// initialisation des 5 paramètres au coeur
core.p0 = tCent; // température du coeur
core.p1 = dCent; // densité
core.p2 = rayonCoeur; // rayon (en rayon solaire)
core.p3 = luminositeCoeur; // Luminosité (en luminosité solaire)
core.p4 = masseCoeur; // masse du coeur (en masse solaire)
iteration = 0;
do {
iteration = iteration +1; // compteur d'itérations
core = integre1pas(core); // on intègre
deviation = Math.abs((core.p4/fitMass)-1.0); // écart à la masse au point de jonction
core.p0 = Math.max(core.p0,tEffective); // pour toujours avoir T>0
core.p1 = Math.max(core.p1,1.0e-12); // pour toujours avoir rho>0
if(plotControl == 1) {
if(plotMasseOuRayon) {
xx = core.p2/radius; // pour des tracés en fonction du rayon
}
else {
xx = core.p4/mass; // pour des tracés en fonction de la masse
}
if((xx >0.001) && (xx < 0.999)) {
vecteurXout.addElement(new Double(xx));
yy = Math.log(core.p0)/Math.log(10);
vecteurTout.addElement(new Double(yy));
yy = Math.log(core.p1)/Math.log(10);
vecteurDout.addElement(new Double(yy));
if(plotMasseOuRayon) {
yy = Math.log(core.p4)/Math.log(10);
}
else {
yy = Math.log(core.p2)/Math.log(10);
}
vecteurRout.addElement(new Double(yy));
yy = Math.log(core.p3)/Math.log(10);
vecteurLout.addElement(new Double(yy));
index = (int) Math.floor(xx*keepSize); // indice de tableau
T_array[index] = core.p0;
D_array[index] = core.p1;
R_array[index] = core.p2;
L_array[index] = core.p3;
M_array[index] = core.p4;
K_array[index] = kOpacite;
E_array[index] = puissance;
P_array[index] = pression;
}
}
}
while((deviation > fitPrecision) && (iteration < maxiteration)); // fin de la boucle do...while
// le point de jonction n'est pas atteint après maxIteration itérations
if( deviation > fitPrecision) {
System.out.println("dehors:trop loin de la jonction. Déviation=" + deviation);
}
return core;
}
/** Intégration de la structure stellaire depuis la surface de l'étoile.
Entrée : température effective et luminosité
Sortie : les 5 valeurs au point de jonction (T, densité, rayon, luminosité, masse)
*/
public Vecteur5 dedans(double tEff, double lum, int error){
Vecteur5 surf = new Vecteur5();
// Constantes(s) locale(s)
// -----------------------
final int maxiteration = 1000; // nombre maximum d'itérations
// Variable(s) locale(s)
// ---------------------
int iteration; // compteur d'itérations
double deviation; // écart de la masse au dernier point à la masse au point de jonction
double xx, yy;
int index;
radius = Math.sqrt(lum)/Math.pow((tEff/SOLEIL_T_SURFACE),2); // rayon de la photosphère (en rayon solaire)
// initialisation des 5 paramètres à la surface
surf.p0 = tEff/twoTo14; // temp de surface plus petite que Teffect
surf.p1 = 1.0e-12; // densité très faible
surf.p2 = radius; // rayon (en rayon solaire)
surf.p3 = lum; // Luminosité (en luminosité solaire)
surf.p4 = mass; // masse du coeur (en masse solaire)
iteration = 0;
do {
iteration = iteration +1; // compteur d'itérations
surf = integre1pas(surf); // on intègre
deviation = Math.abs((surf.p4/fitMass)-1.0); // écart à la masse au point de jonction
surf.p3 = Math.max(surf.p3,1.0e-10); // pour toujours avoir L>0
if(plotControl == 1) {
if(plotMasseOuRayon) {
xx = surf.p2/radius; // pour des tracés en fonction du rayon
}
else {
xx = surf.p4/mass; // pour des tracés en fonction de la masse
}
if((xx > 0.001) && (xx < 0.999)) {
vecteurXin.addElement(new Double(xx));
yy = Math.log(surf.p0)/Math.log(10);
vecteurTin.addElement(new Double(yy));
yy = Math.log(surf.p1)/Math.log(10);
vecteurDin.addElement(new Double(yy));
if(plotMasseOuRayon) {
yy = Math.log(surf.p4)/Math.log(10);
}
else {
yy = Math.log(surf.p2)/Math.log(10);
}
vecteurRin.addElement(new Double(yy));
yy = Math.log(surf.p3)/Math.log(10);
vecteurLin.addElement(new Double(yy));
index = (int) Math.floor(xx*keepSize); // indice de tableau
T_array[index] = surf.p0;
D_array[index] = surf.p1;
R_array[index] = surf.p2;
L_array[index] = surf.p3;
M_array[index] = surf.p4;
K_array[index] = kOpacite;
E_array[index] = puissance;
P_array[index] = pression;
}
}
}
while((deviation > fitPrecision) && (iteration < maxiteration)); // fin de la boucle do...while
// le point de jonction n'est pas atteint après maxIteration itérations
if( deviation > fitPrecision) {
System.out.println("dedans:trop loin de la jonction. Déviation=" + deviation);
}
return surf;
}
/** Correction des conditions aux limites pour l'itération suivante.
Entrée : anciennes conditions aux limites et valeurs des dérivées
Sortie : valeur corrigée pour Teff, densit" et température au coeur, luminosité
Cette méthode utilise la technique de BP pour trouver les corrections.
*/
public void correctionValeur(double deriv[][]){
// Constantes(s) locale(s)
// -----------------------
final int numEqs=4; // nombre d'équations
// Variable(s) locale(s)
// ---------------------
int n=numEqs; // à des fins de comparaison vs. Golub/Van Loan
int mu; // position du pivot d'ordre k
int pivk; // position du pivot au pas k
double buffer;
double factor; // facteur de réduction pour les corrections, ou bien correction max
double maxmat; // valeur max. de la sous-matrice pour determiner mu
// résolution d'un système de 4 équations pour T,rho,r,L: sommation sur j:
// deriv(i,j)*corrections(j)+deriv(i,5)=0
// LU decomposition (Golub/Van Loan p.112) n=NumEqs
mu=0;
for( int k=0 ; k <= n-2 ; k++ ) {
maxmat=0;
for( int i=k ; i <= n-1 ; i++ ) {
for ( int j=0 ; j <= k ; j++ ) {
if (Math.abs(deriv[i][j]) >= maxmat) {
maxmat = Math.abs(deriv[i][j]);
mu=i-k;
}
}
}
pivk = mu + k;
// permutation des Lignes k/pivk
if (pivk != k) {
for ( int i=k ; i <= n ; i++ ) {
buffer = deriv[k][i];
deriv[k][i] = deriv[pivk][i];
deriv[pivk][i] = buffer;
}
}
if (deriv[k][k] != 0.) {
for ( int i=k+1 ; i <= n-1 ; i++ ) {
deriv[i][k] = deriv[i][k]/deriv[k][k];
}
for ( int i=k+1 ; i <= n-1 ; i++ ) {
for ( int j=k+1 ; j <= n ; j++ ) {
deriv[i][j] = deriv[i][j]-deriv[i][k]*deriv[k][j];
}
}
}
else {
System.out.println("correctionValeur: division par 0 dans calcul de deriv[i][j]");
}
}
// on introduit le second membre, qui est -delta(:)
for ( int i=0 ; i <= n-1 ; i++ ) {
corrections[i] = -deriv[i][n];
}
// Back Substitution (Golub/Van Loan p.89)
if ( deriv[n-1][n-1] != 0. ) {
corrections[n-1] = corrections[n-1]/deriv[n-1][n-1];
}
else {
System.out.println("correctionValeur: division par 0 dans calcul de corrections[n-1]");
}
for (int i=n-2 ; i >= 0 ; i-- ) {
double dot_product = 0;
for (int j=i+1 ; j<= n-1 ; j++) {
dot_product=dot_product + deriv[i][j]*corrections[j];
}
if ( deriv[i][i] != 0. ) {
corrections[i] = (corrections[i]-dot_product)/deriv[i][i];
}
else {
System.out.println("correctionValeur: division par 0 dans calcul de corrections[i]");
}
}
// les corrections inconnues del(z(i)) ont été trouvées
// vérification de la taille des corrections
double maxinter = 0;
factor=1;
for (int i=0 ; i <= numEqs-1 ; i++) {
factor = Math.abs(corrections[i]/maxCorrect[i]);
factor = Math.max(factor,maxinter);
maxinter = factor;
factor = Math.max(1.,factor);
}
for (int i=0 ; i <= numEqs-1 ; i++) {
corrections[i] = corrections[i]/factor;
}
// on utilise les del(z(i)) pour corriger les conditions aux limites
tEffective = tEffective*(1+corrections[0]);
luminosite = luminosite*(1+corrections[1]);
tCentral = tCentral*(1+corrections[2]);
dCentral = dCentral*(1+corrections[3]);
}
}
class Vecteur5 {
double p0;
double p1;
double p2;
double p3;
double p4;
}