PRO RANDOM_FIX, MSTAR, DPC, A_FIX, M_FIX, NGEN, OUTGEN, FULL=FULL, $ R_STAR=R_STAR, L_STAR=L_STAR, T_STAR=T_STAR, A_STAR=A_STAR, MG_STAR=MG_STAR, $ BAND=BAND, EFIX=EFIX, OFIX=OFIX, IFIX=IFIX, ALB_FIX=ALB_FIX, $ OBS_DATE=OBS_DATE, JDSTART=JDSTART, TIMESPAN=TIMESPAN ;+ ; NAME: RANDOM_FIX ; PURPOSE: ; EXPLANATION: given the values of mass and period of a ; planet, mass and distance of its host star, this code generates a ; fixed number of orbital parameters and the corresponding physical ; properties ; CALLING SEQUENCE: RANDOM_FIX, MSTAR, DPC, P_FIX, M_FIX, OBS_DATE, $ ; NGEN, OUTGEN ; INPUTS: ; MSTAR - Stellar mass (Msun) ; DPC - distance of the star (pc) ; A_FIX - fixed value of the semi-major axis (AU) of the planets to ; be generated ; M_FIX - fixed value of the mass (Mjup) of the planets to ; be generated ; NGEN - number of planets to be generated ; OPTIONAL INPUT KEYWORDS: ; EFIX - Fixed value of the eccentricity ; OFIX - Fixed value of Omega_node (deg), ; IFIX - Fixed value of the inclination (deg), ; T0FIX - Fixed value of the time of periastron passa_star ; OBS_DATE - obs. date (JD-2400000.0D) ; JDSTART - Starting date of the simulated observations, if ; obs_real is not set ; TIMESPAN - duration of the survey,if obs_real is not set ; FULL - to be set to 1 if the full evaluation of the physical ; caracteristics of the planets is required. In this ; case, the following keywords MUST be provided too: ; R_STAR - Stellar radius (Rsun) ; L_STAR - Bolometric luminosity of the star - cgs units ; A_STAR - age of the star (Myrs) ; MG_STAR - absolute magnitude of the star in the chosen band ; BAND - spectral range for the observations (strarr) ; T_STAR - Stellar Temperature (K) ; ; ; OUTPUTS: ; OUTGEN - .idl file containing all the characteristics of the ; generated planets: ; mass - planet mass, as given as input ; p - planet period, as given as input ; a - semi-major axis, obtained using the 3rd Kepler's law ; e - eccentricity ; irad - orbit inclination (rad) ; omega_peri_rad - long of periastron (rad) ; omega_node_rad - long of node (rad) ; frac_ill - illuminated fraction of ; the planet suface, ; evaluated with projection.pro ; rho_arcsec - projected separation in arcsec ; x_arcsec - x position on the projected orbit ; y_arcsec - y position on the projected orbit ; pla_rad - planet radius (Rjup) ; pla_teff - planet effective temperature (K) ; pla_cnt - planet/star contrast ; kappa_rv - planet RV signature ; ampl_astrom - planet astrometric signature ; RESTRICTIONS: ; ; PROCEDURES CALLED: ; PROJECTION(), REPLICATE(), CONTR_RAD(), INTRINSIC_TEFF(), ; GAS_RADIUS(), ROCKY_RADIUS() ; ; REVISION HISTORY: ; Written M. Bonavita, March, 2009 ; Modified ;- ;#####################################################################; rsun=6.95508E8 ; m Raggio del Sole rjup=7.13E7 ; m Raggio di Giove rear=6.378E6 ; m Raggio della Terra rje=rjup/rear ; Raggio di Giove in unita' di raggi terrestri rsj=rsun/rjup ; Raggio del sole in unita' di raggi gioviani au=1.496E11 ;m Unità Astronomica sigma=(5.67*1.e-5) ;costante di Stefan Boltzman in unità cgs L_sun=3.82*1.e33 ;luminosità solare in unità cgs h_const = 6.63e-34 ; Planck's constant J*s c_light = 3.00e8 ; speed of light m/s k_const = 1.38e-23 ; Boltzmann's constant J/K lambda_c=11.3*1e-6 hc_lk=(h_const*c_light)/(lambda_c*k_const) ;#####################################################################; a=replicate(a_fix,ngen) p=(365.25*a^(3.D/2.D))/(Mstar(0)^(1.D/2.D)) mass=replicate(m_fix,ngen) ;p=replicate(p_fix, ngen) ;a=((p/365.25D)^(2.D/3.D))*mstar^(1.D/3.D) ;;;;;;;;;;;;, eccentricity ;;;;;;;;;;;;;;;;;;;;;;;;;;; e=0.667*RANDOMU(seed,ngen, /double) if keyword_set(efix) then e=replicate(efix(0), ngen) if e[0] gt 1 then e=0.667*RANDOMU(seed,ngen, /double) ;;;;;;;;;;;;, inclination ;;;;;;;;;;;;;;;;;;;;;;;;;;; cosi=2.0*RANDOMU(seed,ngen, /double)-1.0 irad=ACOS(cosi) ideg=irad*!radeg if keyword_set(ifix) then ideg=replicate(ifix(0), ngen) & irad=ideg/!radeg if ideg[0] gt 360 then begin cosi=2.0*RANDOMU(seed,ngen, /double)-1.0 irad=ACOS(cosi) ideg=irad*!radeg end ;;;;;;;;;;;;;;;;;;;;; omega (long of periastron) ;;;;;;;;;;;;;;;;;;;;;;;;;; omega_peri=RANDOMU(seed,ngen, /double)*360.0 omega_peri_rad=omega_peri/!radeg ;;;;;;;;;;;;, Time of periastron pasage ;;;;;;;;;;;;;;;;;;;;;;;;;;; t0=fltarr(ngen) FOR i=0l, ngen-1 DO t0(i)=RANDOMU(seed,1,/double)*p(i) ;days ;;;;;;;;;;;;;;;;;;;;; Omega (long of node) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; omega_node=RANDOMU(seed,ngen, /double)*360.0 if keyword_set(ofix) then omega_node=replicate(ofix(0), ngen) if omega_node[0] gt 360 then omega_node=RANDOMU(seed,ngen, /double)*360.0 omega_node_rad=omega_node/!radeg ;;;;;;;;;;;;;;;;;;;;; Obs. time (if not provided) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; IF keyword_set(obs_date) eq 0 then obstime=RANDOMU(seed, ngen, /double)*timespan*365.25+jdstart IF keyword_set(obs_date) then obstime=replicate(obs_date, ngen) ;;;;;; end of random generation of orbital elements ;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; CALCULATION OF THE PARAMETER OF THE PROJECTED ORBIT ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; rho=DBLARR(ngen) rho_arcsec=DBLARR(ngen) frac_ill=DBLARR(ngen) project=PROJECTION(p, a, e, irad, omega_peri_rad, omega_node_rad, t0, obstime, ngen, dpc) rho_arcsec=project.rho/dpc[0] x_arcsec=project.x2/dpc[0] y_arcsec=project.y2/dpc[0] frac_ill=project.frac_ill IF keyword_set(FULL) eq 0 THEN BEGIN save, ngen, mass, a, p, e, rho_arcsec, x_arcsec,y_arcsec, omega_peri_rad, irad, omega_node_rad, t0, obstime, filename=outgen goto, skip END ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; CALCULATION OF THE PHYSICAL PROPERTIES OF THE PLANETS ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Geometric Albedo pla_alb=0.2+(RANDOMU(seed,ngen, /double)*0.5) IF keyword_set(alb_fix) then pla_alb=replicate(alb_fix[0],ngen) ; Temperature ; Effective temperature (Baraffe et al. 2003) Temp_cool=intrinsic_teff(mass, a_star*1.e-3, band) ;Equivalent temperature (see Sudarsky et al. 2001) Temp_rflx=((1.-pla_alb)*L_Star/(16.*!pi*sigma*((a*(au*1e+2))^2.)))^(1/4.) ; Planet Radius if a_star*1.e-3 ge 1.5 then kind='COLD' if a_star*1.e-3 lt 1.5 and a_star*1.e-3 gt 0.5 then kind='WARM' if a_star*1.e-3 le 0.5 then kind='HOT' fm=0.3 ;ice/rock or rock/iron fraction for rocky planets rtype=round(randomu(seed, 1)) ;type of rocky planets if rtype eq 0 then type='ROCKY' ;the planet is composed by rock and iron if rtype eq 1 then type='ICE' ;the planet is composed by Ice and rock if mass[0] gt 40./318. then pla_rad=replicate(gas_radius(mass[0], 10/318., a[0], kind=kind), ngen) ;jupiter like planets from Fortney et al., 2007 if mass[0] gt 10./318. and mass[0] le 40./318. then pla_rad=replicate(10^(0.0086+0.3755*ALOG10(mass[0])), ngen) ;neptune like planets: solar system mass-radius relation if mass[0] le 10./318. then pla_rad=replicate((rocky_radius(mass[0]*318., fm, type=type))/rje, ngen) ;rocky planets from Fortney et al., 2007 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; CALCULATION OF THE OBSERVABLE QUANTITIES ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Planet/Star contrast IF band NE 'nband' THEN BEGIN pla_cnt=fltarr(ngen) Pla_Teff=((Temp_rflx^4.)+(Temp_cool^4.))^(1./4.) FOR jj=0, ngen-1 DO BEGIN out=contr_rad(mg_star, mstar,a_star*1.e-3,band, mass[0], a[0], frac_ill[jj], pla_rad[0], pla_alb[jj]) pla_cnt[jj]=sqrt((out.contrast_ref^2.)+(out.contrast_int^2.)) IF kind EQ 'COLD' THEN pla_cnt[jj]=out.contrast_ref & Pla_Teff[jj]=Temp_rflx[jj] IF kind EQ 'HOT' THEN pla_cnt[jj]=out.contrast_int & Pla_Teff[jj]=Temp_cool[jj] ;print, frac_ill[jj], pla_alb[jj], pla_cnt[jj];, out.contrast_ref, out.contrast_int, sqrt((out.contrast_ref^2.)+(out.contrast_int^2.)) END ;FOR jj=0, ngen-1 END ;if band ne 'nband' if band eq 'nband' then begin LBB_pl=exp(hc_lk/pla_teff)-1. LBB_st=exp(hc_lk/T_Star)-1. contr_Nband_int=(LBB_st/LBB_pl)*((pla_rad*rjup)/(R_Star*rsun))^2. contr_Nband_ref=1/out.contrast_ref contrast_N=sqrt(contr_Nband_int^2+contr_Nband_ref^2) end ; Radial velocity amplitude p_yr=p/365.25 kappa_rv=28.4*p_yr^(-0.33)*mass*sin(irad)/((mstar^0.667)*SQRT(1-e^2)) ; astrometric signature ampl_astrom=0.960*(a/5.0)*10.0/dpc*mass/mstar ;mas ; transit probability and depth TR_pb=(a*COS(irad))/R_Star TR_dp=pla_rad/R_Star ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; save, ngen, mass, a, p, e, rho_arcsec, x_arcsec,y_arcsec, omega_peri_rad, irad, omega_node_rad, t0, obstime, frac_ill, pla_rad, pla_teff, pla_cnt, kappa_rv, ampl_astrom, TR_pb, TR_dp, pla_alb, filename=outgen skip: end ;PROGRAM