// //++++++++++++++++ //.IDENTIFICATION $Id: dms_simu.java,v 2.0 2002/01/09 17:06:22 arenou Exp $ //.TYPE class //.LANGUAGE java //.AUTHOR F.Arenou //.ENVIRONMENT UNIX and MacOS above 8.1 //.KEYWORDS orbits,Gaia //.PURPOSE class for the generation of a random multiple system //.COMMENTS This class allows to create a random system: given a star //.COMMENTS with some absolute magnitude, the star will become single or //.COMMENTS double, etc, with some probability. If not single, the masses //.COMMENTS and orbital parameters of its companions are chosen at random. //---------------- // import java.util.* ; ///////////////////////////////////////////////////////////////////////////// public class dms_simu extends dms_orbit { ///////////////////////////////////////////////////////////////////////////// protected double [] absmag, vmicol ; // absolute mag, color private double [] tab_param = new double[6] ; // orbital parameters // //++++++++++++++++ //.NAME dms_simu //.PURPOSE generation of a double (or multiple) system //.INPUT absolute magnitude of primary (double) //.INPUT (optional) V-I color of primary (double) //.INPUT (optional) Mass (M_sun) of primary (double) //---------------- // public dms_simu(double MV, double VmI, double Mass) { super() ; // explicit empty construction init_simu(MV,VmI,Mass) ; } // computes roughly the mass assuming a main-sequence star public dms_simu(double MV, double VmI) { super() ; // explicit empty construction init_simu(MV,VmI,Mag2Mass(MV)) ; } // computes roughly mass and V-I assuming a main-sequence star public dms_simu(double MV) { super() ; // explicit empty construction init_simu(MV,BmV2VmI(MV2BmV(MV)),Mag2Mass(MV)) ; } // //++++++++++++++++ //.NAME getNumberOfCompanions //.OUTPUT returns the number of companions (int) //.COMMENTS The number of components is getNumberOfCompanions()+1 //---------------- // public int getNumberOfCompanions() { return number_companions ; } // //++++++++++++++++ //.NAME getAbsMag //.INPUT number (>= 0) of this companion (int) //.OUTPUT returns the absolute magnitude of the component (double) //---------------- // public double getAbsMag(int which_comp) { int comp_nb = init_with_bound(which_comp,0,number_companions, "Companion number") ; return absmag[comp_nb] ; } // //++++++++++++++++ //.NAME getVmICol //.INPUT number (>= 0) of this companion (int) //.OUTPUT returns the V-I color of the component (double) //---------------- // public double getVmICol(int which_comp) { int comp_nb = init_with_bound(which_comp,0,number_companions, "Companion number") ; return vmicol[comp_nb] ; } // // this is the true initialisation, whatever the constructor is // private void init_simu(double MV, double VmI, double star_mass) { number_companions = numb_companions() ; // at random absmag = new double[number_companions+1] ; absmag[0] = MV ; vmicol = new double[number_companions+1] ; vmicol[0] = VmI ; super.primary(number_companions,star_mass) ; // creates the primary if (number_companions == 0) return ; // poor star alone // now get random values of orbital parameters for all orbparam[0][0]=0 ; // is used for a test on period for (int comp_nb=1; comp_nb <= number_companions; comp_nb++) { int this_comp = comp_nb ; double M2, P ; // mass and period if (rand_numb(0.,1.) < 0.1) { // make about 5% planet M2 = mass_planet() ; P = planet_period() ; } else { do { M2 = mass_ratio()*star_mass ;} while (M2 < 0.01) ; // > 10 Jup P = period() ; if ((M2 < 0.08)&&(P < 3000)) { // make short-P BD rare do { P = period() ; } while (P < 3) ; } } if (P < orbparam[comp_nb-1][0]) { // need to swap components this_comp = comp_nb-1 ; super.companion(comp_nb,super.getMass(this_comp), super.getCompanionInfo(this_comp)) ; absmag[comp_nb] = absmag[this_comp] ; vmicol[comp_nb] = vmicol[this_comp] ; } // fill orbital parameters array tab_param[0] = P ; tab_param[1] = rand_numb(0.,tab_param[0]) ; // random T tab_param[2] = eccentricity() ; // random e tab_param[3] = rand_numb(0.,360.) ; // random omega2 double cosi = rand_numb(0.,1.) ; // random i tab_param[4] = Math.atan2(Math.sqrt(1.-cosi*cosi),cosi)/DegToRad ; tab_param[5] = rand_numb(0.,360.) ; // random Omega // now creates the companion with these parameters super.companion(this_comp,M2,tab_param) ; double MV2 = Mass2Mag(M2) ; // absolute magnitude if (MV > MV2) MV2=MV ; // should never happen absmag[this_comp] = MV2 ; // abs. mag and vmicol[this_comp] = BmV2VmI(MV2BmV(MV2)) ; // color of companion } } // // generation of the number of stellar companions // Duquesnoy&Mayor 91, uncorrected for incompleteness // but adding 6% of planets+BD // private int numb_companions() { int rand = rand_numb(0,106) ; if (rand < 51) return 0 ; else if (rand < 51+40 +3) return 1 ; else if (rand < 51+40+7 +5) return 2 ; else return 3 ; // 2 +1% } // // random generation of the period (days) // Duquesnoy&Mayor 91 // private double period() { double log10p=rand_gaus(4.8,2.3,-1,10) ; // bell-shaped return(Math.pow(10.,log10p)) ; } // // random generation of the planet period (days) // private double planet_period() { double log10p=rand_numb(0.48,5) ; // constant? return(Math.pow(10.,log10p)) ; } // // random generation of the mass ratio // we do not discuss whether or not there is a peak at q=1 // private double mass_ratio() { double q=rand_numb(0.,1.) ; if (q < 0.1) { // increasing up to 0.1 if (q/0.1 < rand_numb(0.,1.)) q=rand_numb(0.,1.) ; } // else constant return(q) ; } // // random generation of the planet mass // in 1/M between about 1 and 10 M_Jup // private double mass_planet() { double log10M=rand_numb(-3.,-2.) ; return(Math.pow(10.,log10M)) ; } // // random generation of the excentricity // from Duquesnoy&Mayor 91 // private double eccentricity() { if (tab_param[0] < 10) return 0 ; // circularization below 10 days else if (tab_param[0] < 1000) return rand_gaus(0.33,0.2,0.,1.) ; // bell-shaped else return Math.sqrt(rand_numb(0.,1.)) ; // f(e)=2e above 1000 days } // // mass-luminosity relation Henry & Mac Carthy 93 // mass (solar mass) --> absolute magnitude // private double HMMass2Mag(double mass) { double lum, logM ; logM=log10(mass) ; if (mass > 2) lum=3.78-5.9*logM ; // ca, c'est moi qui l'a ajoute else if (mass > 0.5) lum=19.77-0.5*Math.sqrt(1628.66*logM+852.49) ; else if (mass > 0.18) lum=(1.4217-logM)/0.1681 ; else if (mass > 0.08) lum=22.36-0.5*Math.sqrt(760.89*logM+925.31) ; else lum = 17.59 ; return(lum) ; } // // mass-luminosity relation // mass (solar mass) --> absolute magnitude // probably wrong at the bottom of MS and surely wrong // for BD and planets, but this will not be used // so we just need to avoid an overflow // private double Mass2Mag(double mass) { return(4.77-2.5*log10(Mass2Lum(mass))) ; } private double Mass2Lum(double mass) { if (mass <= 0) { return(1.e-30) ; } else if (mass <= 0.08) { return(mpow(2063,mass,9)) ; } else { // Tout et al., 1996, MNRAS 281, 257 (3% precision on luminosity) return (mpow(0.3970417,mass,5.5)+ mpow(8.527626,mass,11))/( mpow(0.00025546,mass,0)+ mpow(1,mass,3)+ mpow(5.432889,mass,5.2)+ mpow(5.563579,mass,7)+ mpow(0.7886606,mass,8)+ mpow(0.00586685,mass,9.5)) ; } } // // mass-luminosity relation Henry & Mac Carthy 93 // absolute magnitude --> mass (solar mass) // private double Mag2Mass(double mag) { double logM ; if (mag < 1.45) logM=(3.78-mag)/5.9 ; // ca, c'est moi qui l'ai ajoute else if (mag < 10.25) logM=0.002456*mag*mag-0.09711*mag+0.4365 ; else if (mag < 12.89) logM=-0.1681*mag+1.4217 ; else { if (mag > 17.59) mag=17.59 ; logM=0.005257*mag*mag-0.2351*mag+1.4124 ; } return(Math.pow(10.,logM)) ; } // // gives B-V as a fct of MV (Henry AJ 512) // !!! wrong for BD/EP // private double MV2BmV(double x) { double [] T = {-1.242359,0.555204,-0.038184,0.000992} ; if (x < -10) x = -10 ; if (x > 20) x = 20 ; return(T[0]+x*(T[1]+x*(T[2]+x*(T[3])))) ; } // // gives (V-I)c as a fct of (B-V)J // !!! wrong for BD/EP // private double BmV2VmI(double x) { double [] T = {0.02129,1.3119,-0.14209,-1.333,1.5637,-0.3924} ; if (x < -0.3) x = -0.3 ; return(T[0]+x*(T[1]+x*(T[2]+x*(T[3]+x*(T[4]+x*(T[5])))))) ; } // // gives a random number between [min and max[ // private double rand_numb(double min, double max) { return min+(max-min)*Math.random() ; } private int rand_numb(int min, int max) { return min+(int)((max-min)*Math.random()) ; } // // poors man gaussian, truncated outside of [min,max] // private double rand_gaus(double moy, double sig, double min, double max) { double res ; do { res=-6 ; for (int i=1; i<=12; i++) res += Math.random() ; res *= sig ; res += moy ; } while ((res < min)||(res > max)) ; return res ; } // // decimal log // private double log10(double x) { return Math.log(x)/ln10 ; } }