// //++++++++++++++++ //.IDENTIFICATION $Id: dms_orbit.java,v 2.0 2002/01/09 17:04:28 arenou Exp $ //.TYPE class //.LANGUAGE java //.AUTHOR F.Arenou //.ENVIRONMENT UNIX and MacOS above 8.1 //.KEYWORDS orbits,Gaia //.PURPOSE Computes the orbits for all components of a system. //.COMMENTS Once the orbital parameters are known, the position offset //.COMMENTS (in A.U.), radial vel., flux, LLT of each component may be known //.COMMENTS Assumes *no* interaction between companions. //.COMMENTS for the moment, we assume that the motion of a companion n //.COMMENTS is around the barycenter of the n-1 inner companions //.COMMENTS and the effects on the primary are just added //---------------- // import java.text.* ; ///////////////////////////////////////////////////////////////////////////// public class dms_orbit { ///////////////////////////////////////////////////////////////////////////// // constants during the life of this system; here [] means for all components protected int number_companions ; // # of companions (apart of primary) protected double [] mass, radius ; // mass (solar), radius (A.U) protected double [][] orbparam ; // orbital parameters of all but first protected int [] PrimaryNb ; // this is the number of primary protected double [] MassRatio ; // with respect to primary protected double [] SecondAxis ; // second. semi-maj. axis of all but 1st protected double [] SecondAmpli ; // semi-amplitude of secondary (km/s) protected double LastEpoch ; // last time the position was computed // variable arrays, depending of observation epoch protected double [] ksi, eta ; // position on tangent plane (A.U.) protected double [] vra, ltt ; // radial velocity and LTT effect protected double [] rad, flx ; // radius vector, relative flux (eclipses) // gives min and max limits for some parameters protected int[] limit_companion={0,4} ; // # of comp protected double[] limit_mass={(float)1.e-6,150} ; // the mass protected double[] limit_radius={(float)1.e-6,20} ; // the radius protected double[] limit_perd={(float)0.1,(float)1.e+10} ; // period, days protected double[] limit_exce={(float)0.,(float)0.999999} ; // excentricity // shortcuts in order to avoid computations protected double [] sinw, cosw, sinOmg, cosOmg, sinI, cosI, efac ; protected double Pi=3.141592653589793, DegToRad=Pi/180.0, ln10=Math.log(10.) ; protected double AUMkm = 149.597870691 ; // 1 U.A. in Gm protected double RSunMkm = 0.696265, RSunAU = RSunMkm/AUMkm ; // sun radius // //++++++++++++++++ //.NAME dms_orbit //.PURPOSE constructs all the components of the multiple system //.INPUT number of companions to the primary (int) //.INPUT masses (in solar mass) of all components (double [numberofcomp+1]) //.INPUT orbital parameters of all components (double [numberofcomp+1][6]) //.INPUT (optional) reference primary (default=0) (int [numberofcomp+1]) //.INPUT (optional) radius (in solar radius) of all components (double [numberofcomp+1]) //.COMMENTS orbital param = P (days), T (days), e, omega2 (deg), i (deg), Omega (deg) //.COMMENTS orbital_parameters[0] is not used //---------------- // public dms_orbit(int nb_companions, double [] init_masses, double [][] init_orb_param, int [] reference, double [] init_radii) { primary(nb_companions, init_masses[0], init_radii[0]) ; // init primary for (int comp=1; comp<=nb_companions; comp++) // other comps companion(comp, init_masses[comp], init_orb_param[comp], reference[comp], init_radii[comp]) ; } public dms_orbit(int nb_companions, double [] init_masses, double [][] init_orb_param, int [] reference) { primary(nb_companions, init_masses[0]) ; // init primary for (int comp=1; comp<=nb_companions; comp++) // other comps companion(comp, init_masses[comp], init_orb_param[comp], reference[comp]) ; } public dms_orbit(int nb_companions, double [] init_masses, double [][] init_orb_param) { primary(nb_companions, init_masses[0]) ; // init primary for (int comp=1; comp<=nb_companions; comp++) // other comps companion(comp, init_masses[comp], init_orb_param[comp]) ; } //empty construction should be followed by a call to "primary" and "companion" public dms_orbit() { number_companions=-1 ; } // take no risk // //++++++++++++++++ //.NAME primary //.PURPOSE constructs the primary of the multiple system and allocate data //.INPUT number of companions to the primary (int) //.INPUT mass (in solar mass) of primary (double) //.INPUT (optional) radius (in solar radius) (double) //---------------- // public void primary(int companions, double mass_primary, double radius_primary) { make_primary(companions, mass_primary, radius_primary*RSunAU) ; } public void primary(int companions, double mass_primary) { make_primary(companions, mass_primary, Mass2Radius(mass_primary)) ; } // //++++++++++++++++ //.NAME companion //.PURPOSE creates one of the companions of the multiple system //.INPUT number (> 0) of this companion (int) //.INPUT mass (in solar mass) of this companion (double) //.INPUT orbital parameters (see below) (double [6]) //.INPUT (optional) number of reference primary (int) =0 by default //.INPUT (optional) radius (in solar radius) (double) //.COMMENTS P (days), T (days), e, omega2 (deg), i (deg), Omega (deg) //---------------- // public void companion(int init_number, double init_mass, double[] init_orb_param, int reference, double init_radius) { make_companion(init_number,init_mass,init_orb_param,reference,init_radius*RSunAU) ; } public void companion(int init_number, double init_mass, double[] init_orb_param, int reference) { make_companion(init_number,init_mass,init_orb_param,reference,Mass2Radius(init_mass)) ; } public void companion(int init_number, double init_mass, double[] init_orb_param) { make_companion(init_number,init_mass,init_orb_param,0,Mass2Radius(init_mass)) ; } // //++++++++++++++++ //.NAME printCompanionInfo //.PURPOSE gives information about the orbital parameters of the component //.INPUT number (> 0) of this companion (int) //---------------- // public void printCompanionInfo(int which_comp) { DecimalFormat fmt = new DecimalFormat("######.######"); String[][] info = { {"P (days)", "period"}, {"T (days)", "time of periastron passage"}, {"e (-)", "eccentricity"}, {"o2 (deg)", "argument of periastron"}, {"i (deg)", "inclination"}, {"O (deg)", "position angle of the node"}, } ; int comp_nb = init_with_bound(which_comp,1,number_companions, "Companion number") ; System.out.println("\t" + "semi-major axis" + " " + "a (A.U.)" + " = " + fmt.format(SecondAxis[comp_nb]) + " (reflex=" + fmt.format(SecondAxis[comp_nb]*MassRatio[comp_nb]) + ")") ; System.out.println("\t" + "semi-amplitude" + " " + "K (km/s)" + " = " + fmt.format(SecondAmpli[comp_nb]) + " (reflex=" + fmt.format(SecondAmpli[comp_nb]*MassRatio[comp_nb]) + ")") ; for (int i=0; i<6; i++) { double temp = orbparam[comp_nb][i] ; if (i > 2) temp /= DegToRad ; // angles in degrees System.out.println("\t" + info[i][1] + " " + info[i][0] + " = " + temp) ; } } // //++++++++++++++++ //.NAME getCompanionInfo //.INPUT number (> 0) of this companion (int) //.OUTPUT the orbital parameters of the component (double [6]) //.COMMENTS the angles are in radian, not degrees //---------------- // public double [] getCompanionInfo(int which_comp) { int comp_nb = init_with_bound(which_comp,0,number_companions, "Companion number") ; return orbparam[comp_nb] ; } // //++++++++++++++++ //.NAME Orbital_Position //.PURPOSE computes the position of all components at some date //.INPUT date (days) (double) //.COMMENTS the date is for instance with respect to HJD 2440000, but //.COMMENTS doesnt matter, provided it is consistent with T periastron date //---------------- // public void OrbitalPosition(double time) { // initialization of primary for (int i=0; i<=number_companions ; i++) { FillOutput(i,0.,0.,0.,0.,0.) ; flx[i] = 0 ; } LastEpoch = time ; if (number_companions != 0) Do_Orbit(time) ; // else do nothing } // //++++++++++++++++ //.NAME getMass //.INPUT number (>= 0) of this companion (int) //.OUTPUT the mass of the component (double) //---------------- // public double getMass(int which_comp) { int comp_nb = init_with_bound(which_comp,0,number_companions, "Companion number") ; return mass[comp_nb] ; } // //++++++++++++++++ //.NAME getRadius //.INPUT number (>= 0) of this companion (int) //.OUTPUT returns the radius (R_sun) of the component (double) //---------------- // public double getRadius(int which_comp) { int comp_nb = init_with_bound(which_comp,0,number_companions, "Companion number") ; return radius[comp_nb]/RSunAU ; } // //++++++++++++++++ //.NAME isEclipsing //.INPUT number (> 0) of this companion (int) //.OUTPUT returns true if eclipsing binary (boolean) //.COMMENTS assuming null eccentricity //---------------- // public boolean isEclipsing(int which_comp) { int comp_nb = init_with_bound(which_comp,1,number_companions, "Companion number") ; double a = SecondAxis[comp_nb]*(1+MassRatio[comp_nb]) ; // relative semi-major axis return (a*Math.abs(cosI[comp_nb]) < (radius[PrimaryNb[comp_nb]]+radius[comp_nb])) ; } // //++++++++++++++++ //.NAME get_ksi //.PURPOSE gives the current x position on the tangent plane of a component //.INPUT number (>= 0) of this companion (int) //.OUTPUT x offset in A.U. wrt the barycenter of the system (double) //---------------- // public double getKsi(int star_nb) { int comp_nb = init_with_bound(star_nb,0,number_companions, "Companion number") ; return ksi[comp_nb] ; } // //++++++++++++++++ //.NAME get_eta //.PURPOSE gives the current y position of a component on the tangent plane //.INPUT number (>= 0) of this companion (int) //.OUTPUT y offset in A.U. wrt the barycenter of the system (double) //---------------- // public double getEta(int star_nb) { int comp_nb = init_with_bound(star_nb,0,number_companions, "Companion number") ; return eta[comp_nb] ; } // //++++++++++++++++ //.NAME get_rho //.PURPOSE gives the separation between the component and the primary //.INPUT number (>= 0) of this companion (int) //.OUTPUT separation (A.U.) between component and primary (double) //---------------- // public double getRho(int star_nb) { int comp_nb = init_with_bound(star_nb,0,number_companions, "Companion number") ; double dx = ksi[comp_nb]-ksi[PrimaryNb[comp_nb]] ; double dy = eta[comp_nb]-eta[PrimaryNb[comp_nb]] ; return(Math.sqrt(dx*dx+dy*dy)) ; } // //++++++++++++++++ //.NAME get_rad //.PURPOSE gives the distance offset along line of sight of a component //.INPUT number (>= 0) of this companion (int) //.OUTPUT distance offset in A.U. wrt the barycenter (double) //---------------- // public double getRad(int star_nb) { int comp_nb = init_with_bound(star_nb,0,number_companions, "Companion number") ; return rad[comp_nb] ; } // //++++++++++++++++ //.NAME get_vra //.PURPOSE gives the radial velocity offset of a component //.INPUT number (>= 0) of this companion (int) //.OUTPUT radial velocity offset in km/s wrt the barycenter (double) //---------------- // public double getVra(int star_nb) { int comp_nb = init_with_bound(star_nb,0,number_companions, "Companion number") ; return vra[comp_nb] ; } // //++++++++++++++++ //.NAME get_ltt //.PURPOSE gives the light-time travel effect of a component //.INPUT number (>= 0) of this companion (int) //.OUTPUT LTT offset in days wrt the barycenter (double) //---------------- // public double getLtt(int star_nb) { int comp_nb = init_with_bound(star_nb,0,number_companions, "Companion number") ; return ltt[comp_nb] ; } // //++++++++++++++++ //.NAME get_flx //.PURPOSE gives the flux factor of a component //.INPUT number (>= 0) of this companion (int) //.OUTPUT flux factor (1 outside eclipse, [0-1[ during ecl.) (double) //---------------- // public double getFlx(int star_nb) { int comp_nb = init_with_bound(star_nb,0,number_companions, "Companion number") ; return flx[comp_nb] ; } // // True initialisation of primary // private void make_primary(int init_companions, double init_mass_primary, double init_radius_primary) { number_companions = init_with_bound(init_companions,limit_companion[0], limit_companion[1],"Number of companions") ; GetRoom() ; // now allocating enough ram for all data // init parameters of primary mass[0] = init_with_bound(init_mass_primary,limit_mass[0],limit_mass[1], "Primary mass (in M_sun)") ; // in solar mass MassRatio[0]=1 ; SecondAxis[0]=0 ; SecondAmpli[0]=0 ; // for completeness radius[0] = init_with_bound(init_radius_primary,limit_radius[0],limit_radius[1], "Primary radius (in R_sun)") ; // in solar radius for (int i=0; i<=5; i++) orbparam[0][i] = 0. ; // no orbital param } // // True initialisation of one companion // private void make_companion(int init_number, double init_mass, double[] init_orb_param, int init_ref, double init_radius) { int comp_nb = init_with_bound(init_number,1,number_companions, "Companion number") ; int ref_nb = init_with_bound(init_ref,0,number_companions-1, "Reference number") ; // number of the primary for this star if (ref_nb > 0) { System.out.println("Other reference component than 0 "+ "not yet implemented") ; ref_nb=0 ; } PrimaryNb[comp_nb] = ref_nb ; mass[comp_nb] = init_with_bound(init_mass,limit_mass[0],limit_mass[1], "Companion mass (in M_sun)") ; radius[comp_nb] = init_with_bound(init_radius,limit_radius[0],limit_radius[1], "Companion radius (in R_sun)") ; // in solar radius make_orbparam(comp_nb, init_orb_param) ; // the orbital parameters make_shortcuts(comp_nb) ; // a few shortcuts to avoid trigonometry at each epoch make_axis(comp_nb) ; // semi-major axis and semi-amplitude } // // initialisation of the orbital parameters // private void make_orbparam(int comp_nb, double[] init_orb_param) { double [] fP = orbparam[comp_nb] ; // a shortcut for the parameter table fP[0] = init_with_bound(init_orb_param[0],limit_perd[0],limit_perd[1], "Orbital period (in days)") ; fP[1] = init_orb_param[1] ; // time of periastron fP[2] = init_with_bound(init_orb_param[2],limit_exce[0],limit_exce[1], "Excentricity") ; fP[3] = (init_orb_param[3] % 360)*DegToRad ; // omega_2 in radians fP[4] = (init_orb_param[4] % 180)*DegToRad ; // i fP[5] = (init_orb_param[5] % 360)*DegToRad ; // Omega } // // initialisation of annex parameters // private void make_shortcuts(int comp_nb) { double [] fP = orbparam[comp_nb] ; // a shortcut for the parameter table sinw[comp_nb] = Math.sin(fP[3]); cosw[comp_nb] = Math.cos(fP[3]); sinI[comp_nb] = Math.sin(fP[4]); cosI[comp_nb] = Math.cos(fP[4]); sinOmg[comp_nb] = Math.sin(fP[5]) ; cosOmg[comp_nb] = Math.cos(fP[5]); efac[comp_nb] = Math.sqrt((1.+fP[2])/(1.-fP[2])) ; } // // initialisation of semi-major axis // private void make_axis(int comp_nb) { double [] fP = orbparam[comp_nb] ; // a shortcut for the parameter table double ref_mass = 0 ; for (int i=0; i 0; comp_nb--) { Orbit_Component(time, comp_nb) ; } // check if eclipse, with possible LTT effect for (int comp_nb=1; comp_nb<=number_companions ; comp_nb++) { if (isEclipsing(comp_nb)) Eclipse(time, comp_nb) ; } // convert relative eclipsed surface to flux loss for (int i=0; i<=number_companions ; i++) { // ? several simultaneous eclipses ;-) flx[i] = 1-flx[i] ; if (flx[i] < 0) flx[i]=0 ; } } // // this computes the orbit position for one companion // private void Orbit_Component(double time, int comp_nb) { double [] pos = new double[2] ; // ksi, eta double [] fP = orbparam[comp_nb] ; // a shorcut for the parameter table double E = Kepler(fP[2], 2.*Pi*(time-fP[1])/fP[0]) ; double v = 2.*Math.atan2(efac[comp_nb]*Math.sin(E/2.),Math.cos(E/2.)) ; double Q=v+fP[3]; double oneminusecosE = 1-fP[2]*Math.cos(E) ; // first, for the component itself OrbPos(SecondAxis[comp_nb]*oneminusecosE, Q, comp_nb, pos) ; ksi[comp_nb] += pos[0] ; eta[comp_nb] += pos[1] ; rad[comp_nb] += SecondAxis[comp_nb]*oneminusecosE*Math.sin(Q)*sinI[comp_nb] ; vra[comp_nb] += VrPos(SecondAmpli[comp_nb],Math.cos(Q),fP[2],cosw[comp_nb]) ; ltt[comp_nb] += LttPos(SecondAxis[comp_nb],v,Math.sin(Q),fP[2],comp_nb) ; // then just add the reflex motion, for the inner bodies Q -= Pi ; double PrimAxis = MassRatio[comp_nb]*SecondAxis[comp_nb] ; double PrimAmpli = MassRatio[comp_nb]*SecondAmpli[comp_nb] ; OrbPos(PrimAxis*oneminusecosE, Q, comp_nb, pos) ; double rad0 = PrimAxis*oneminusecosE*Math.sin(Q)*sinI[comp_nb] ; double vra0 = VrPos(PrimAmpli,Math.cos(Q),fP[2],cosw[comp_nb]) ; double ltt0 = LttPos(PrimAxis, v, Math.sin(Q), fP[2], comp_nb) ; for (int primary=0; primary < comp_nb; primary++) { ksi[primary] += pos[0] ; eta[primary] += pos[1] ; rad[primary] += rad0 ; vra[primary] += vra0 ; ltt[primary] += ltt0 ; } } // // coordinates of current orbit points // private void OrbPos(double r, double Q, int comp_nb, double [] coo ) { double sinQ = Math.sin(Q); double cosQ = Math.cos(Q); ; coo[0] = r*(cosQ*sinOmg[comp_nb] + sinQ*cosOmg[comp_nb]*cosI[comp_nb]) ; coo[1] = r*(cosQ*cosOmg[comp_nb] - sinQ*sinOmg[comp_nb]*cosI[comp_nb]) ; } // // radial velocity for current position // private double VrPos(double K, double CosQ, double e, double Cosw) { return ((CosQ+e*Cosw)*K) ; } // // light time (days) for current position // private double LttPos(double a, double v, double sinQ, double e, int comp_nb) { double cosv = Math.cos(v); ; return (AUMkm*a*sinI[comp_nb]/25900)* ((1-e*e)*sinQ/(1+e*cosv)+e*sinw[comp_nb]) ; } // // light dimming due to eclipse at current position // from P. North, Cours de Goutelas, 2000 // on output, flx is the relative surface, not the dimming (=1-relative_surface) // private void Eclipse(double time, int comp_nb) { double delta ; // apparent sep in A.U // take LTT effect of a third component into account if (number_companions > 1) { double [][] bck = new double[number_companions+1][5] ; // backup for (int i=0; i<=comp_nb; i++) { bck[i][0]=ksi[i]; bck[i][1]=eta[i]; bck[i][2]=rad[i]; bck[i][3]=vra[i]; bck[i][4]=ltt[i] ; } FillOutput(0,0.,0.,0.) ; FillOutput(comp_nb,0.,0.,0.) ; Orbit_Component(time-bck[0][4], comp_nb) ; delta = getRho(comp_nb) ; for (int i=0; i<=comp_nb; i++) { FillOutput(i,bck[i][0],bck[i][1],bck[i][2],bck[i][3],bck[i][4]) ; // restore } } else { delta = getRho(comp_nb) ; } if (delta > radius[PrimaryNb[comp_nb]]+radius[comp_nb]) { // no eclipse flx[comp_nb] = 0 ; } else { double a = SecondAxis[comp_nb]*(1+MassRatio[comp_nb]) ; delta /= a ; // d/(a1+a2) => no unit double sintheta2 = (delta*delta-cosI[comp_nb]*cosI[comp_nb])/ (sinI[comp_nb]*sinI[comp_nb]) ; double theta = Math.atan2(Math.sqrt(sintheta2),Math.sqrt(1-sintheta2)) ; double r1 = radius[0]/a, r2 = radius[comp_nb]/a, k = r2/r1 ; double sinrho1=0, cosrho1=0, sinrho2=0.5, cosrho2=0 ; for (int i=0; i<30; i++) { double lastsinrho2 = sinrho2 ; sinrho1=k*sinrho2 ; cosrho1 = Math.sqrt(1-sinrho1*sinrho1) ; cosrho2=(delta/r1-cosrho1)/k ; if (cosrho2 > 1) cosrho2=1 ; if (cosrho2 < -1) cosrho2=-1 ; sinrho2 = Math.sqrt(1-cosrho2*cosrho2) ; if (Math.abs(lastsinrho2-sinrho2) < 0.001) break ; } double rho1 = Math.atan2(sinrho1,cosrho1) ; double rho2 = Math.atan2(sinrho2,cosrho2) ; double alpha = ((rho1-cosrho1*sinrho1)+k*k*(rho2-cosrho2*sinrho2))/Pi ; if (rad[comp_nb] < rad[0]) { flx[0] += alpha ; flx[comp_nb] = 0 ;} // companion in front else { flx[comp_nb] = alpha/(k*k) ;} // host in front } return ; } // // Fill output parameters // private void FillOutput(int comp_nb, double ksi0, double eta0, double rad0, double vra0, double ltt0) { ksi[comp_nb]=ksi0; eta[comp_nb]=eta0; rad[comp_nb]=rad0; vra[comp_nb]=vra0; ltt[comp_nb]=ltt0; } private void FillOutput(int comp_nb, double ksi0, double eta0, double rad0) { ksi[comp_nb]=ksi0; eta[comp_nb]=eta0; rad[comp_nb]=rad0; } // // Solve Kepler's Equation // private double Kepler(double e, double M) { double E = M + e*Math.sin(M) + e*e*Math.sin(2*M)/2 ;// eccentric anomaly for(int j=0; j<10; j++) E += (M-E+e*Math.sin(E))/(1-e*Math.cos(E)); return E; } // // mass-radius relation // mass (solar mass) --> radius (A.U.) // private double Mass2Radius(double mass) { double InSolarRadius ; if (mass <= 0) { InSolarRadius = 0.; // for BD and planets, little interest. !!! We do not yet account // on larger radius for hot jupiters } else if (mass <= 0.045) { InSolarRadius = 0.1 ; } else if (mass <= 0.08) { InSolarRadius = 0.135-0.875*mass ; } else { // Tout et al., 1996, MNRAS 281, 257 (1.2% precision on radius) InSolarRadius = (mpow(1.715359,mass,2.5)+ mpow(6.597788,mass,6.5)+ mpow(10.08855,mass,11)+ mpow(1.012495,mass,19)+ mpow(0.07490166,mass,19.5))/( mpow(0.01077422,mass,0)+ mpow(3.082234,mass,2)+ mpow(17.84778,mass,8.5)+ mpow(1,mass,18.5)+ mpow(0.00022582,mass,19.5)); } return InSolarRadius*RSunAU ; } // // allocates memory for the whole orbital system // OK, something should be done for systems with 1 component only, // although for Gaia this only means several dozen of gigabytes wasted ;-) // private void GetRoom() { PrimaryNb = new int[number_companions+1] ; // who is the primary, mass = new double[number_companions+1] ; // mass in M_sun, radius = new double[number_companions+1] ; // radius in A.U. orbparam = new double[number_companions+1][6] ; // orb. param (0 unused) MassRatio = new double[number_companions+1] ; // M2/M1 SecondAxis = new double[number_companions+1] ; // a2 (0 is unused) SecondAmpli = new double[number_companions+1] ; // K2 (0 is unused) ksi = new double[number_companions+1] ; // position in x, eta = new double[number_companions+1] ; // position in y, rad = new double[number_companions+1] ; // radial velocity, vra = new double[number_companions+1] ; // radial velocity, ltt = new double[number_companions+1] ; // light, flx = new double[number_companions+1] ; // flux, sinOmg = new double[number_companions+1] ; // +a few shortcuts... cosOmg = new double[number_companions+1] ; sinI = new double[number_companions+1] ; cosI = new double[number_companions+1] ; sinw = new double[number_companions+1] ; cosw = new double[number_companions+1] ; efac = new double[number_companions+1] ; } // // coefficient times mass at some power // protected double mpow(double coef, double mass, double exponent) { return(coef*Math.exp(exponent*Math.log(mass))) ; } // // initialises a value with a check of min and max value // protected float init_with_bound(float init, float min, float max, String message) { float value ; if ((init < min) || (init > max)) { if (init < min) value=min ; else value=max ; System.out.println(message + " out of range [ " + min + " - " + max + " ], value " + value + " assumed.") ; } else value=init ; return value ; } protected double init_with_bound(double init, double min, double max, String message) { double value ; if ((init < min) || (init > max)) { if (init < min) value=min ; else value=max ; System.out.println(message + " out of range [ " + min + " - " + max + " ], value " + value + " assumed.") ; } else value=init ; return value ; } protected int init_with_bound(int init, int min, int max, String message) { int value ; if ((init < min) || (init > max)) { if (init < min) value=min ; else value=max ; System.out.println(message + " out of range [ " + min + " - " + max + " ], value " + value + " assumed.") ; } else value=init ; return value ; } }