Program driveInterL !C. Mordasini, University of Bern, Switzerland, 2019 !Simple driver program for the subroutine to get the interpolated luminosity !of low-mass planets with H/He envelopes as a function of time. implicit none double precision :: Mc,Me,t,Lumi,L1,L2,L3,L4,L5,L6 Integer ::i !Get core mass, envelope mass, and age from user input Write(*,*)'Enter core mass in Mearth (1-20)' READ(*,*)Mc Write(*,*)'Enter envelope mass in Mearth (1-20)' READ(*,*)Me Write(*,*)'Enter planet age in yrs (1e8-1e10)' READ(*,*)t !Call the interpolation subroutine to get the luminosity in LJupiter call InterL(Mc,Me,t,Lumi) !Print the result WRITE(*,*)'Luminosity in intrinsic LJupiter',Lumi WRITE(*,*)'Luminosity in intrinsic LEarth ',Lumi/0.0000896588 READ(*,*) !Make some tables in fort.14 and fort.15 files !1) temoral evolution for the Mcore and Me entered by the user Do i=1,990 t=1e8+10**(REAL(i)/100D0) call InterL(Mc,Me,t,Lumi) !WRITE(*,*)t,Mc,Me,Lumi WRITE(14,*)t,Mc,Me,Lumi end do !2) temporal evoluton for the Mcore entered by the user !and a choice of 6 different envelope masses Do i=1,990 t=1e8+10**(REAL(i)/100D0) Me=0 call InterL(Mc,Me,t,L1) Me=1e-4 call InterL(Mc,Me,t,L2) Me=1e-2 call InterL(Mc,Me,t,L3) Me=1e-1 call InterL(Mc,Me,t,L4) Me=1e0 call InterL(Mc,Me,t,L5) Me=1e1 call InterL(Mc,Me,t,L6) !WRITE(*,*)t,L1,L2,L3,L4,L5,L6 WRITE(15,*)t,L1,L2,L3,L4,L5,L6 end do end Program driveInterL subroutine InterL(Mc,Me,t,Lumi) !C. Mordasini, University of Bern, Switzerland, 2018 ! !Subroutine to estimate the intrinsic luminosity of distant low-mass planets with H/He envelopes as !a function of time by interpolation in fitting functions (C. Mordasini, A&A, 2019) derived from !evolutionary calculations with the COMPLETO model. The core has an Earth-like composition with !1/3 iron and 2/3 silicates (no ices) !Input !Mc (double precision): planet core mass in units of Mearth. Allowed domain: 1-20 !Me (double precision): planet H/He envelope mass in units of Mearth. Allowed domain: 0-20 !t (double precision): planet age in years. Allowed domain: 1e8 to 1e10 ! !Output !Lumi (double precision): planet intrinsic luminosity in units of Jupiter's intrinsic luminosity today LJ (LJ=8.67e-10 Lsun) implicit none Double Precision,INTENT(IN) ::Mc,Me,t Double Precision,INTENT(OUT) ::Lumi Integer ::i,il,iu Double Precision ::tg,Ll,Lu Logical,PARAMETER ::debug=.FALSE. INTEGER,PARAMETER ::Nttab=14,Npar=5 Double Precision,Dimension(Nttab) ::ttab=(/0.1D0, 0.3D0, 0.5D0, 0.8D0, 1D0, 2D0, 3D0, 4D0, 5D0, 6D0, 7D0, 8D0, 9D0, 10D0/) Double precision,PARAMETER ::Mcmax=20D0,Memax=20D0 Double Precision,Dimension(Nttab,Npar) ::fc !The fitting parameters fc(1,:)=(/0.00252906,-0.00200238,0.00104408,0.0586485,0.000967878/) fc(2,:)=(/0.00121324,-0.000533601,0.000360703,0.0214114,0.000368533/) fc(3,:)=(/0.000707416,-0.000394131,0.000212475,0.0138138,0.000189456/) fc(4,:)=(/0.000423376,-0.000187283,0.000125872,0.00887292,0.000117141/) fc(5,:)=(/0.000352187,-0.00014148,9.94382e-05,0.00718831,9.20563e-05/) fc(6,:)=(/0.000175775,-4.07832e-05,4.5853e-05,0.00357941,5.52851e-05/) fc(7,:)=(/0.00011412,-2.09944e-05,2.91169e-05,0.00232693,4.00546e-05/) fc(8,:)=(/8.81462e-05,-2.07673e-05,2.12932e-05,0.00171412,2.90984e-05/) fc(9,:)=(/6.91819e-05,-1.90159e-05,1.62128e-05,0.00134355,2.30387e-05/) fc(10,:)=(/5.49615e-05,-1.6862e-05,1.29045e-05,0.00109019,1.96163e-05/) fc(11,:)=(/4.5032e-05,-1.51951e-05,1.05948e-05,0.00091005,1.70934e-05/) fc(12,:)=(/3.80363e-05,-1.40113e-05,8.93639e-06,0.00077687,1.50107e-05/) fc(13,:)=(/3.30102e-05,-1.31146e-05,7.69121e-06,0.000675243,1.32482e-05/) fc(14,:)=(/2.92937e-05,-1.24023e-05,6.73922e-06,0.000595191,1.17809e-05/) !convert age in Gyr tg=t/1D9 !Check if in tabulated domain in time and Mc and Me IF(tgttab(Nttab))THEN WRITE(*,*)'Age higher then maximal tabulated value',tg,ttab(Nttab) WRITE(*,*)'STOPPING' STOP END IF IF(Mc>Mcmax)THEN WRITE(*,*)'Mcore higher then max fitted value',Mc,Mcmax WRITE(*,*)'STOPPING' READ(*,*) END IF IF(Me>Memax)THEN WRITE(*,*)'Menve higher then max fitted value',Me,Memax WRITE(*,*)'STOPPING' READ(*,*) END IF !get time index Do i=1,Nttab-1 IF(ttab(i)<=tg.AND.ttab(i+1)>=tg)THEN il=i iu=i+1 EXIT END IF end do IF(DEBUG)WRITE(*,*)il,iu,ttab(il),ttab(iu) !do fitting Ll=fc(il,1)+fc(il,2)*Mc+fc(il,3)*Mc**2D0+fc(il,4)*Me+fc(il,5)*Me**2D0 Lu=fc(iu,1)+fc(iu,2)*Mc+fc(iu,3)*Mc**2D0+fc(iu,4)*Me+fc(iu,5)*Me**2D0 IF(DEBUG)WRITE(*,*)Ll,Lu !interpolate linearly in log(L(log(t))) CALL LinIntPolCI(log10(tg),log10(Ll),log10(ttab(il)),log10(Lu),log10(ttab(iu)),Lumi) !convert back from log(L) to L for output lumi=10**lumi Return End subroutine InterL SUBROUTINE LinIntPolCI(x,y1,x1,y2,x2,y) !makes linear interpolation for f(x) (Dim n) at the point x, if (x1,f(x1)=y1) and (x2,f(x2)=y2) !are given (dimension y1, y2 also n of course) using Lagrange's formula IMPLICIT NONE DOUBLE PRECISION ::y1,y2,y DOUBLE PRECISION ::x,x1,x2 LOGICAL,PARAMETER ::debug=.FALSE. IF(debug)THEN WRITE(*,*) 'Linintpol x,y1,y2,x1,x2,y',x,y1,y2,x1,x2,y END IF y=((x-x2)/(x1-x2))*y1+((x-x1)/(x2-x1))*y2 IF(debug)WRITE(*,*) 'y after',y END SUBROUTINE LinIntPolCI