2. 1. planetesimal dust and gas accretion test, under alpha phase
4. alpha experimental semi-scientific only, simple assumptions and code
5. 1. python 3 code
6. 1. 03.05.2024 0000.06.07.14

import os import sys import getopt import math import numpy as np from scipy.stats import norm import pandas as pd

import matplotlib.pyplot as plt

1. 1. some constants

pi=math.pi degrad=pi/180

year_in_seconds=31556926 seconds_in_year=year_in_seconds megayear1=year_in_seconds*1e6

G=6.6743015e-11 kb=1.380649e-23 sb= 5.670374419e-8 planck=6.62607015e-34 rydberg=10973731.6

avogadro=6.02214076e23 mproton=1.67262192e-27 dalton=1.660539066605e-27 ## amu amu=dalton

au=149597870700 ## meters

teff_sun=5772 tsun=teff_sun msun=1.988471e30 lsun=3.828e26 rsun=6.957e8

mearth=5.9722e24 rearth=6.3781e6 rooearth=5517

mjupiter=1.89813e27 mjupiter_mearth=317.828

def accrete_mass_from_dustgas_1(mstar1,snowline1, incdustcoeff1, incgascoeff1, bhill1, time1, timestep, taugas, tausolids, aas1, eccs1, ms1, mrocks1, mices1, mgases1, diska0, dust1, gas1, ices1, rocks1): ## accrete dust and gas mstar_me1=330000*mstar1 bee=2.4*1 ## 2.4 3 3.4 ... aas2=np.copy(aas1) eccs2=np.copy(eccs1) ms2=np.copy(ms1) mgases2=np.copy(mgases1) mices2=np.copy(mices1) mrocks2=np.copy(mrocks1) dust2=np.copy(dust1) gas2=np.copy(gas1) ices2=np.copy(ices1) rocks2=np.copy(rocks1)

miso1=10*mstar1 ming1=5 # 10 minimum gas accretion mass gapmass1=30 #print("dt",timestep) len1=len(aas1)

for n in range(0, len1):

a1=aas1[n] m1=ms1[n]

mpgas1=mgases1[n] mpice1=mices1[n] mprock1=mrocks1[n]

e1=eccs1[n] ## location in disk array

idx1 = (np.abs(diska0 - a1)).argmin() ## note dust disk, not planetesimals mdust1=dust1[idx1] rock1=rocks1[idx1] ice1=ices1[idx1] mgas1=gas1[idx1]

#iceprop1=ice1/(rock1+ice1) #rockprop1=rock1/(rock1+ice1) #iceprop1=ice1/mdust1 #rockprop1=rock1/mdust1 rockprop1=1 iceprop1=0 if(a1>snowline1): rockprop1=0.5 iceprop1=0.5

aas2[n]=a1

rhill1=math.pow((m1/(mstar_me1*3)),1/3) omega1=math.sqrt(mstar1*np.power(a1, -3/2)) r1=math.pow(m1, 0.27) vesc1=math.sqrt(m1/r1) vorb1=math.sqrt(mstar_me1/a1) vhill1=math.sqrt(m1/rhill1) #sucklen1=rhill1*3.5 sucklen1=rhill1*bhill1 ## original 3.5 sucklen0=rhill1*bhill1 suckdist_est0=mdust1*sucklen0 sdensedust0=mdust1 ## approx sdenseplanet0=m1/sucklen0 ## sigma*mass accretion_runaway_switch_1=1 ## check this if( (2*sdenseplanet0*m1) > (sdensedust0*(mdust1*sucklen0)) ): accretion_runaway_switch_1=0

runaway1=math.pow(m1, 4/3) oligark1=math.pow(m1, 2/3)

#dustaddfactor1=0.01 #dustaddfactor1=0.00000066 dustaddfactor1=incdustcoeff1*1.7e-6 ## 1.25 e-6 ## 6 e-7 maddfactor1=dustaddfactor1*timestep*math.pow(mstar1,1) #dustaddmk1=mdust1*sucklen1*omega1*maddfactor1 ## dust accretion dustaddmk1=mdust1*sucklen1*maddfactor1 ## dust accretion

if (m1<miso1): #print("under miso", m1) ms2[n]=m1+dustaddmk1 mrocks2[n]=mprock1+dustaddmk1*rockprop1 mices2[n]=mpice1+dustaddmk1*iceprop1 else: #print("over miso", m1) ms2[n]=m1+dustaddmk1*1e-7 ## mininal add on isolation mass! mrocks2[n]=mprock1+dustaddmk1*rockprop1*1e-7 mices2[n]=mpice1+dustaddmk1*iceprop1*1e-7

#gasaddfactor1=0.20 #gasaddfactor1=incgascoeff1*0.028 gasaddfactor1=incgascoeff1*5e-3 ## 0.65 0.1

if(m1<gapmass1): ## no gap formed gasu1=gasaddfactor1*timestep*math.pow(rhill1, 3)*math.pow(mstar1, 1) else: ## gap gasu1=gasaddfactor1*timestep*math.pow(rhill1, 2)*math.pow(mstar1, 1)

if(m1>5000): mummoo=1 #print("m1 ming1: ", m1, ming1)

if(m1>ming1): ## gap formation gapk1=0.1 ## gap depth from hat estim if(m1>gapmass1): gapk1=0.015 if(m1>(gapmass1*10)): gapk1=0.001 if(m1>(gapmass1*30)): gapk1=0.0001

gasu1=gasu1*gapk1

if(m1>2000): gasu1=0

ms2[n]=m1+gasu1 mgases2[n]=mpgas1+gasu1 #print("GAS",mgases2[n], gasu1)

return(aas2, eccs2, ms2, mrocks2, mices2, mgases2, dust2, gas2, ices2, rocks2)

def calc_macc_1(mstar1, taudisk1, time1): age_myr1=time1*1e-6 macc1=math.pow(10, -8.24)*math.pow(mstar1, 2.4) ## alternative agexp -1.4 ... -2.68 , pow ((age/myr), -2) macc2=macc1*np.exp(-age_myr1/taudisk1) return(macc2)

def calculate_sigmas_from_macc_0(mstar1, macc1, aas1, alpha1): disk_q=-1/2 ## temperature exponnet #https://staff.fnwi.uva.nl/c.dominik/Teaching/SPF/07-viscous_accretion.pdf year_in_seconds=31556926 msun=1.988471e30 mproton=1.67262192e-27 au=149597870700 ## meters G=6.6743015e-11 kb=1.380649e-23 muu1=2.34 ## from hat molecular h/he + macc2=macc1*(msun/year_in_seconds) mstar2=mstar1*msun aas2=np.copy(aas1)*au #alpha1=0.01/(macc1/1e-7) ## check this tdisk1=np.copy(aas2) ## CHECK THIS tempdisk1=650*np.power((aas2/au), disk_q) ## must obtain better! sigmagases1=((muu1*mproton*math.sqrt(G*mstar2))/(3*math.pi*kb))*macc2/(alpha1*tempdisk1*np.power(aas2, 3/2)) return(sigmagases1)

def generate_povray(name1,aas1, rs1, planettypes1, drawskala1): filename1="aplanets1.pov" starname1=name1

names1=["a","b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m", "n","o"]

file1 = open(filename1, "w")

lines1=["#include \"colors.inc\" \n" "#include \"textures.inc\" \n" "#include \"stones.inc\" \n" "#include \"functions.inc\" \n" "\nlight_source {<1000, -1000, 2000> color rgb <1,1,1> *2 }\n" "\ncamera {location <25, -100, 200> look_at <25, 0, 0> right 5 up 1 angle 15}\n" "\n object { text { ttf \"timrom.ttf\" \""+starname1+"\" 0.1, 0 } scale 3 "+ "translate <20,3,0>\n "+"texture {pigment { color rgb 1} } }" "\nsphere { <0, 0, 0>, 0.33 }"]

file1.writelines(lines1)

xk1=6*drawskala1 xofset1=0

rk1=1 ## radius render koeff

len1=len(aas1) for n in range(0,len1): mm=n type1=planettypes1[n] #px1=math.log10(aas1[n]*100)-2*ak1 a1=aas1[n] #la1=math.log10(aas1[n]*100) #px1=la1*xk1+xofset1 px1=a1*xk1+xofset1

#pr1=math.log10(rs1[n])*rk1 pr1=math.sqrt(rs1[n])*rk1 #if(type1==1): pr1=pr1*3 #pr1=1 print(n, aas1[n], rs1[n], planettypes1[n]) #print(n,px1,pr1) if(type1==0): ## base lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment {color rgb <1,0.9,0.7>*0.5 }} }\n"] if(type1==1): ## stone lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { function { f_wrinkles(x,y,z) } color_map { [0 color rgb <1,0.8,0.5>] [1 color rgb <1,1,1>*0.5]} }} }\n"] if(type1==2): ## ice lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment {color rgb <1,1,1>*0.85 }} }\n"] if(type1==3): ## neptune lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { function { y+f_wrinkles(y,1,2) } sine_wave frequency 1 scale 4 warp { turbulence 0.5 } scale 0.25 color_map { [0 color rgb <0.5,0.9,1>] [1 color rgb <0.7,1,1>*0.7]} }} }\n"] if(type1==4): ## saturn lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { function { y+f_wrinkles(y,1,2) } sine_wave frequency 1 scale 4 warp { turbulence 0.5 } scale 0.25 color_map { [0 color rgb <1,0.9,0.7>] [1 color rgb <0.6,0.5,0.3>*0.9]} }} texture { pigment { gradient y frequency 2 sine_wave color_map { [0 color rgbt 1] [1 color rgbt <1,1,1,1> ] } } } }\n"] if(type1==5): ## jupiter lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { function { y+f_wrinkles(y,1,2) } sine_wave frequency 1 scale 4 warp { turbulence 0.5 } scale 0.25 color_map { [0 color rgb <0.713725, 0.517647, 0.415686>] [1 color rgb <0.984314, 0.992157, 0.945098>]} }} }\n"] if(type1==10): ## stone lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { function { f_wrinkles(x,y,z) } color_map { [0 color rgb <0.505882, 0.317647, 0.215686>] [1 color rgb <1,1,1>*0.5]} }} }\n"] if(type1==11): ## venusian # venusian lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { granite turbulence 0.2 color_map { [0 color rgb <1,0.8,0.5>] [1 color rgb <1,1,1>*0.8 ]} }} texture { pigment { wrinkles color_map { [0 color rgbt 1] [1 color rgbt <1,1,1,0> ] } } } }\n"] if(type1==12): ## terra-like terran lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { agate scale 2 color_map { [0 color rgb <0, 1, 0>] [0.5 color rgb <0, 1, 0>][0.5 color rgb <0, 0, 1>] [1 color rgb <0, 0, 1>] }} } texture { pigment { wrinkles color_map { [0 color rgbt 1] [1 color rgbt <1,1,1,0> ] } } } }\n"] if(type1==14): ## warm terran moist greenhouse moist warm terran lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { agate scale 2 color_map { [0 color rgb <0.3, 1, 0.2>] [0.5 color rgb <0.3, 1, 0.2>][0.5 color rgb <0.2, 0.2, 1>] [1 color rgb <0.2, 0.2, 1>] }} } texture { pigment { wrinkles scale 0.3 color_map { [0 color rgbt 1] [1 color rgbt <1,1,1,0> ] } } } }\n"] if(type1==15): ## near terran near terran lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { agate scale 2 color_map { [0 color rgb <1,0.8,0.5>] [0.5 color rgb <1,0.8,0.5>][0.5 color rgb <0.5, 0.5, 1>] [1 color rgb <0.5, 0.5, 1>] }} } texture { pigment { wrinkles scale 0.3 color_map { [0 color rgbt 1] [1 color rgbt <1,1,1,0> ] } } } }\n"]

if(type1==13): ## martian martian lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { function { f_wrinkles(x,y,z) } color_map {[0 color rgb <5,0.5,0.5>] [0.3 color rgb <1,0.8,0.5>] [0.9 color rgb <1,1,1>*0.5] [1 color rgb <1,1,1>*1] } }} }\n"] if(type1==20): ## venusian waterworld lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { function { f_wrinkles(1,1,1) } scale 0.3 color_map { [0 color rgb <1, 1, 1>] [1 color rgb <1, 1, 1>]} }} }\n"] if(type1==21): ## ocean waterworld lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { wrinkles scale y/3 color_map { [0 color rgb <1,1,1>] [1 color rgb <0,0,1>]} }} }\n"] if(type1==22): ## icy lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { function { f_wrinkles(x,y,z) } color_map { [0 color rgb <0.8,0.8,0.8>] [1 color rgb <1,1,1>*0.5]} }} }\n"] if(type1==32): ## gasdwarf lines1=["\nobject { sphere {0,"+str(pr1)+"} \n "+"translate <"+str(px1) +",0,0>\n "+"texture {pigment { function { y+f_wrinkles(x,y,z)*0.2 } sine_wave color_map { [0 color rgb <0.7,0.7,1>] [1 color rgb <0.7,0.9,1.0>*1]} }} }\n"]

file1.writelines(lines1) lines1="\n object { text { ttf \"timrom.ttf\" \""+names1[mm]+"\" 0.25, 0 } scale 3 "+ "translate <"+str(px1) +",-5,0>\n "+"texture {pigment { color rgb 1} } }" file1.writelines(lines1)

file1.close()

os.system("povray aplanets1.pov -W3500 -H700 -Q11 -A0.3") #quit(-1) return(0)

def calc_bhill(mstar, mplanet, a, b): msun_me=332831 rhill=a*math.pow( (mplanet/(3*mstar*msun_me)) , 1/3) bhill=b*rhill return(bhill)

def calc_mutual_hill_radius(mstar, a1, m1, a2, m2): msun_me=332831 rhill_m=math.pow( ((m1+m2)/(3*mstar*msun_me)) , 1/3)*((a1+a2)/2) return(rhill_m)

def calc_stability_type_of_inner(a1,a2, m1, m2): stability=0

bhk1=2*math.sqrt(3)

bhill21=calc_bhill(1, m2, a2, bhk1) bhill22=calc_bhill(1, m2, a2, bhk1*2.99) bhill23=calc_bhill(1, m2, a2, bhk1*3.10) ## 3.3 alimit21=a2-bhill21 alimit22=a2-bhill22 alimit23=a2-bhill23 if(a1<alimit21): stability=3 if(a1<alimit22): stability=2 if(a1<alimit23): stability=1

alimit21=alimit21 alimit22=alimit22 alimit23=alimit23

#print(bhill21, alimit21 ) ## asteroids outers #print(bhill22, alimit22 ) ## mars sized #print(bhill23, alimit23 ) ## earth sized

return(stability)

def stability_mutual(mstar1, a1,a2, m1, m2): da1=abs(a2-a1)

mutk1=9 ## 9, 18 ?? rhmut1=calc_mutual_hill_radius(mstar1, a1, m1, a2, m2)

#hillratio1=calc_bhill(mstar1, m1, a1, 1)/calc_bhill(mstar1, m2, a2, 1)

#print("hillratio 1 ", hillratio1)

stable1=0 if(da1>(mutk1*rhmut1)): stable1=1

return(stable1)

def calc_stabilities(mstar1, ajup1, ejup1, mjup1,aas1, eccs1, ms1) :

## calculate jupiter- like planet perturbation power to ## inner planets ## simple hill radius approach ## according of solar system

## default assumpticon stable

stabilities1=np.copy(aas1)*0+1 #0: unsiable #1: stable like earth #2: asteroid zone type stability

bhk1=2*math.sqrt(3)

len1=len(aas1)

for n in range(0, len1):

a1=aas1[n] e1=eccs1[n] m1=ms1[n] da1=ajup1-a1 dm1=m1/mjup1 um1=m1/(mstar1*332831) umjup1=mjup1/(mstar1*332831)

if(a1<ajup1): stability1=calc_stability_type_of_inner(a1,ajup1, m1, mjup1) stabilities1[n]=stability1

return(stabilities1)

def approximate_eccentricities_from_masses(aas1, eccs1, ms1): ## from hat eccentricty approximation by mass of planet

eccs2=np.copy(eccs1) eccs2=np.power((ms1/300),1/1)*0.05 eccs3=np.power((1/ms1),0.75)*0.013 eccs3=np.where(eccs3>0.25, 0.25, eccs3)

len1=len(eccs3)

for n in range(0,len1): if (ms1[n]>10.0): eccs3[n]=eccs2[n]

return(eccs3)

def process_inner_planet_system_by_stability(mstar1, aas1, eccs1, ms1, mrocks1, mices1, mgases1): ## stability v 0001

len1=len(aas1)

aas2=np.copy(aas1) eccs2=np.copy(eccs1) ms2=np.copy(ms1) mrocks2=np.copy(mrocks1) mices2=np.copy(mices1) mgases2=np.copy(mgases1)

## experimental code

eccs2=approximate_eccentricities_from_masses(aas2, eccs2, ms2)

## assumes that "jupiter" only perturbs inner planets

#maxdex1=np.argmax(ms1)

## maximum mass planet in planet system is"jupiter"

jdex1=np.where(ms1>10)[0][1]

for n in range(0,(jdex1-1)): aas1[n]=aas1[n]+eccs1[n] for n in range(jdex1, len(aas1)): aas1[n]=aas1[n]-eccs1[n]

jdexes1=np.where(ms1>10)[0]

#print(ajup1, ejup1, mjup1) #print(jdexes1)

stabilities1=np.copy(aas1)*0+1 stabilities2=np.copy(stabilities1)*0+1

for mm in range((len(jdexes1)-1), -1, -1): #print(mm) jdex1=int(jdexes1[mm]) #jdex1=6 #print(" JUP dex ", jdex1) ajup1=aas1[jdex1] ejup1=eccs1[jdex1] mjup1=ms1[jdex1] #print(ajup1, ejup1, mjup1)

stabilities2=calc_stabilities(mstar1, ajup1, ejup1, mjup1,aas1, eccs1, ms1) len2=len(stabilities2) #print(stabilities2)

for n in range(0, len2):

if (stabilities2[n]==0): stabilities1[n]=stabilities2[n]

if (stabilities2[n]==3): stabilities1[n]=stabilities2[n]

if (stabilities2[n]==2): juppi1=0 if(stabilities1[n]==1): juppi=1 if(juppi1==1): stabilities1[n]=stabilities2[n]

jdex1=int(jdexes1[0]) ajup1=aas1[jdex1] ejup1=eccs1[jdex1] mjup1=ms1[jdex1] #stabilities2=calc_stabilities(mstar1, ajup1, ejup1, mjup1,aas1, eccs1, ms1) #len1=len(aas1)

#stabilities1=np.copy(stabilities2)

#print("===") #print(stabilities1) #print(stabilities2)

#quit(-1)

for n in range(0, len1): ms2[n]=ms2[n] mrocks2[n]=mrocks2[n] mices2[n]=mices2[n] mgases2[n]=mgases2[n]

if (ms1[n]<10): if (stabilities1[n]==0): ms2[n]=0 mrocks2[n]=0 mices2[n]=0 mgases2[n]=0 if (stabilities1[n]==1): ms2[n]=ms1[n] mrocks2[n]=mrocks1[n] mices2[n]=mices1[n] mgases2[n]=mgases1[n] if (stabilities1[n]==2): mrel1=ms1[n]/0.1 ms2[n]=ms1[n]/mrel1 mrocks2[n]=mrocks1[n]/mrel1 mices2[n]=mices1[n]/mrel1 mgases2[n]=mgases1[n]/mrel1 ## can nowcontain gas any ??? if (stabilities1[n]==3): mrel1=ms1[n]/1e-4 ms2[n]=ms1[n]/mrel1 mrocks2[n]=mrocks1[n]/mrel1 mices2[n]=mices1[n]/mrel1 mgases2[n]=0

eccs2=approximate_eccentricities_from_masses(aas2, eccs2, ms2)

len1=len(ms2)

zerodexes1=np.where(ms2==0)[0]

print(zerodexes1)

len1=len(aas1) zlen3=len(zerodexes1)

aas3=[] eccs3=[] ms3=[] mrocks3=[] mices3=[] mgases3=[]

for n in range(0, len1):

mass1=ms2[n]

if(mass1==0): mummu1=0 else: aas3.append(aas2[n]) eccs3.append(eccs2[n]) ms3.append(ms2[n]) mrocks3.append(mrocks2[n]) mices3.append(mices2[n]) mgases3.append(mgases2[n])

#quit(-1)

aas3=np.array(aas3) eccs3=np.array(eccs3) ms3=np.array(ms3) mrocks3=np.array(mrocks3) mices3=np.array(mices3) mgases3=np.array(mgases3)

return(aas3, eccs3, ms3, mrocks3, mices3, mgases3)

def create_ring(mstar1, diskin1, diskout1, diskdelta1, gas1au, amu, asigma, f, asnowline1, gasamountk1, solidamountk1): fc_stonefe=0.25 fc_icefe=1-fc_stonefe fc_iceaddk=(fc_icefe-fc_stonefe)/fc_stonefe

num1=int((diskout1-diskin1)/diskdelta1) adisk1=np.linspace(diskin1, diskout1, num1) len1=len(adisk1) lanesinners1=adisk1-diskdelta1/2 lanesouters1=adisk1+diskdelta1/2 lanesareas1=(math.pi*np.power(lanesouters1, 2)- math.pi*np.power(lanesinners1, 2))*au*au

#sigmagases0=gas1au*np.power(aas1, exp1) ## kg m2 sigmagases0=adisk1*0 gasadd0=norm.pdf(adisk1, loc=amu, scale =asigma)*gas1au

sigmagases0=sigmagases0+gasadd0

sigmagases1=sigmagases0*gasamountk1

sigmasolids_teor1=sigmagases0*f*solidamountk1

snowy1=adisk1

#dustadd1=norm.pdf(adisk1, loc=a1, scale =dustecc1)*dustk1 # mu, sigma = 1, 0.05 # mean and standard deviation #noise1 = np.random.normal(mu, sigma, slices1) #plt.plot(s) #gapp1=1-norm.pdf(aas1, loc=5, scale =0.5) #cavity1=norm.cdf(aas1, loc=5, scale =0.5) #plt.plot(aas1, cavity1,'r-', lw=5, alpha=0.6, label='norm pdf') #plt.show()

snowyk1=norm.cdf(adisk1, loc=snowline1, scale=asnowline1/10) sigmarocks1=sigmasolids_teor1*fc_stonefe iceaddmax1=sigmarocks1*fc_iceaddk sigmaices1=snowyk1*iceaddmax1 sigmasolids1=sigmarocks1+sigmaices1

#plt.plot(sigmasolids1) #plt.plot(sigmaices1)

#plt.show() #quit(-1)

#snowy1np.where(aas1<asnowline1,fc_stonefe,1)

sigmasolids1=snowy1*sigmasolids_teor1

sigmarocks1=sigmasolids_teor1*fc_stonefe sigmaices1=sigmasolids1-sigmarocks1

lanesmasses1=lanesareas1*sigmasolids1 disksolidsmass1=np.sum(lanesmasses1)/msun

print(" Solids ME", disksolidsmass1*33300)

return( adisk1, sigmagases1, sigmasolids1, sigmarocks1, sigmaices1)

def create_disk(mstar1, diskin1, diskout1, diskdelta1, gas1au, exp1, f, asnowline1, gasamountk1, solidamountk1): fc_stonefe=0.25 fc_icefe=1-fc_stonefe fc_iceaddk=(fc_icefe-fc_stonefe)/fc_stonefe

num1=int((diskout1-diskin1)/diskdelta1) aas1=np.linspace(diskin1, diskout1, num1) len1=len(aas1) lanesinners1=aas1-diskdelta1/2 lanesouters1=aas1+diskdelta1/2 lanesareas1=(math.pi*np.power(lanesouters1, 2)- math.pi*np.power(lanesinners1, 2))*au*au

sigmagases0=gas1au*np.power(aas1, exp1) ## kg m2

sigmagases1=sigmagases0*gasamountk1

sigmasolids_teor1=sigmagases0*f*solidamountk1

snowy1=aas1

#dustadd1=norm.pdf(adisk1, loc=a1, scale =dustecc1)*dustk1 # mu, sigma = 1, 0.05 # mean and standard deviation #noise1 = np.random.normal(mu, sigma, slices1) #plt.plot(s) #gapp1=1-norm.pdf(aas1, loc=5, scale =0.5) #cavity1=norm.cdf(aas1, loc=5, scale =0.5) #plt.plot(aas1, cavity1,'r-', lw=5, alpha=0.6, label='norm pdf') #plt.show()

snowyk1=norm.cdf(aas1, loc=asnowline1, scale=asnowline1/10) snowy1=np.where(snowyk1==0,0,1)

sigmaices1=sigmasolids_teor1*snowyk1 sigmarocks1=sigmasolids_teor1*fc_stonefe sigmasolids1=sigmarocks1+sigmaices1

#plt.plot(sigmasolids1) #plt.plot(sigmaices1) #plt.plot(sigmarocks1) #plt.show() #quit(-1)

lanesmasses1=lanesareas1*sigmasolids1 disksolidsmass1=np.sum(lanesmasses1)/msun

#print(" Solids ME", disksolidsmass1*33300) #quit(-1)

return( aas1, sigmagases1, sigmasolids1, sigmarocks1, sigmaices1)

def disk_ring_dust_amount(mstar1, as1, ain1, aout1, sigmagases1, sigmasolids1, sigmarocks1, sigmaices1):

lim1=np.where(aas1>ain1)[0][0] lim2=np.where(aas1>aout1)[0][0]

lanelen1=diskdelta1*au ## inner planets ring

testsigmas1=sigmasolids1[lim1:lim2] testaas1=aas1[lim1:lim2] testlanesinners1=testaas1-diskdelta1/2 testlanesouters1=testaas1+diskdelta1/2 testlanesareas1=(math.pi*np.power(testlanesouters1, 2)- math.pi*np.power(testlanesinners1, 2))*au*au testlanesmasses1=testlanesareas1*testsigmas1

testmass1=np.sum(testlanesmasses1)/mearth return(testmass1)

def greenhouse_from_composition(): #magick=0.7062 #v=TEff?? #tvmagick=tv*magick #teff=tvmagick+tvmagick*21*dno*0.099+tvmagick*dhe*0.218+tvmagick/262*mco2+tvmagick/262*mgc*C return(0)

def calc_greenhouse_tau_co2_h2o(co2pressure, h2opressure): #https://worldbuilding.stackexchange.com/questions/233081/how-to-calculate-greenhouse-effect-for-planet-as-well-as-for-the-atmosphere tau=0.025*np.power(co2pressure*113000, 0.53)+0.277*np.power(h2opressure*113000, 0.3) ## from hat tau=tau*0.016 ## try adjust! return(tau)

def ebm_co2_temperature_degc(Sok, alpha, CO2): Q=(1361.5/4)*Sok ## spherical fastrot sb = 5.670374419e-8 dt = 60. * 60. * 24. * 365. c_w = 4.186E3 rho_w = 1E3 H = 70. #c_w = 0.85E3 #rho_w = 2.5E3 #H = 2 C = c_w * rho_w * H T=0

for n in range(0,300): ASR=(1-alpha)*Q l = np.log(CO2/300.) A = -326.400 + 9.16100*l - 3.16400*l**2 + 0.546800*l**3 B = 1.953 - 0.04866*l + 0.01309*l**2 - 0.002577*l**3 OLR = A + B * T T=T + (dt / C) * (ASR-OLR)

return(T-273.15)

def ebm_greenhouse_estimation(Sok, alpha, pCO2, pressure): ## from hat estimation pco2 absolute, pressure #CO2=pCO2/280e-6 CO2=pCO2/1e-6 T0=ebm_co2_temperature_degc(Sok, alpha, 1) T1=ebm_co2_temperature_degc(Sok, alpha, CO2) T2=T0+(T1-T0)*math.sqrt(pressure) return(T2)

def effective_rau(rau, ecc): rau_eff=rau*math.pow((1-ecc*ecc), 1/4) return(rau_eff)

def calculate_relative_gh_from_p_1(press): ## maybe nok no_habitable_pressure_limit=15e-3 ## bar pressures=[0.1,3] inner_edge_habitable=[0.87, 0.77] outer_edge_habitable=[1.02, 1.18] ## must use log temp to interpolate! #log p vs log t pressus=np.array([0.1,1, 3]) temps_1363watts=np.array([250,288, 310] ) greenhouses_1363=temps_1363watts-255 logpressus=np.log(pressus) logtemps=np.log(temps_1363watts) rp=logpressus-math.log(0.1) rtemps=logtemps-math.log(250) relp1=3 relgh1=55./33. logap=math.log(relp1) logagh=math.log(relgh1) relag=logagh/logap pepe1=math.log(press) pepe2=pepe1*relag geehoo_DT=math.exp(pepe2)*33 return(geehoo_DT)

def calculate_relative_gh_from_p_2(tbase, press): ## only s=1363 W/m2 ## based data from

no_habitable_pressure_limit=15e-3 ## bar pressures=[0.1,3] inner_edge_habitable=[0.87, 0.77] outer_edge_habitable=[1.02, 1.18] ## must use log temp to interpolate! #log p vs log t pressus=np.array([0.1,1, 3]) temps_1363watts=np.array([250,288, 310] ) greenhouses_1363=temps_1363watts-255 logpressus=np.log(pressus) logtemps=np.log(temps_1363watts) rp=logpressus-math.log(0.1) rtemps=logtemps-math.log(250) relp1=3 relgh1=310/288 logap=math.log(relp1) logagh=math.log(relgh1) relag=logagh/logap pepe1=math.log(press) pepe2=pepe1*relag geehoo_T=math.exp(pepe2)*tbase return(geehoo_T)

def t_eq_cloudless(S1, a_ground=0.125):

   cs=13.25e-5 ## in universe, tbackground=2.72 K
   emiss=0.94 ## in most cases bare soli, vegeteed 0.985
   sb=5.670374419e-8
   tground=(2/5)*math.pow((((S1+cs)*(1-a_ground))/(emiss*sb)),0.25)    
   return(tground)

def atmosphere_thermal_enhancement_simplistic(pbar, pco2): co2_ppm=pco2*1e6 ate1=math.pow(pbar, 1/2)*33 ## ate forr earth is 33 dT1=-0.3+2.5*math.log(co2_ppm/280) ## but often dT0=-0.2+1.6*math.log(co2_ppm/280) ate=ate1+dT1 return(ate)

def atmosphere_thermal_enhancement(psurf): ## nte, ate or ghe # ratio between real temperature and gloudless ground temperature ## tex earth cluodless albedo is 0.125 ## empirical relation from solar system ate=math.exp((0.233001*math.pow(psurf, 0.0651203))+(0.0015393*math.pow(psurf, 0.385232))) return(ate)

def temperature_potential_temperature(pressure_kpa):

tpot_temperature=math.pow((pressure_kpa/100), 0.285) return(tpot_temperature)

def water_wapor_pressure_ft(T_degc): wp=611*math.exp((17.27*T_degc)/(237.3+T_degc)) return(wp)

def co2_degc_add(part_co2): co2_ppm=part_co2*1e6 ## odten radiative forcing 5.35*math.log(part_co2/280e-6) ## steeper, maybe realistic dT=-0.3+2.5*math.log(co2_ppm/280) ## but often dT0=-0.2+1.6*math.log(co2_ppm/280) return(dT)

def planet_surface_temperature_from_pressure_official(S1, psurf): ate=atmosphere_thermal_enhancement(psurf) ts_planet=25.3966*math.pow((S1+0.0001325),0.25)*ate return(ts_planet)

def co2_taddition_experimental_1(pco2): pco2_orig=280*1e-6 t_orig=13.7 ## now 380 ppm, 14.4 ## or 1.0 ++ -- 300 -> 400 ppm , 1.4 C 280 --> 400 ppm ## 5x co2, t=3 deg ## note co2 only co2 x2 --0.43c but totally ca 1.9 K water wapor etc() ## because in 150 yrs co2 increased only 0.2 K but actually 0.9 K

dt_co2_2=2.1

num_doublings=math.log(co2_amount/pco2_orig) dt=num_doublings*dt_co2_2 t2=t_input+dt print(t_input, dt, t2) return(dt)

def relative_tau_grom_g(P,gee):

## optical depth #tau=(3*k*P)/(2*gee) reltau=P/gee return(reltau)

def ts_from_teq(teq, tau):

ts=math.pow((1+(tau/2)), 1/4) return(ts)

def rel_temperature_diffu_to_earth(pressure, cp,mol,relomega, relfh2o):

cp0=1000 pressure0=1 #fh2o0=0.0155 #relfh2o=hf2o/2 #relomega=

deerat_base=(pressure/pressure0)*(cp/cp0)*math.pow((mol/29), 2)*math.pow(relfh2o,2) deerat_lat=deerat_base*math.pow(relomega,2) deerat_lon=deerat_base*relomega return(deerat_lat)

def rotperiod_planetmass_density(m,r): ## in earth mss ro=m/math.pow(r,3) prottah=0.632*math.pow(m,1/6)*math.pow((1/ro), 1/3) ## prottah=((2*pi*r)/13.1)*sqrt(m/mjup) ## prottah=pow(mp, 0.27) return(prottah)

def greenhouse_pressu_lab1(pressu1): ## venus, earth, mars actual_temps=[737,288,208] greehouse_tplus=[507,33,8] ### mars 6 K ? atmpressure=[92,1,0.0065] ppres1=92/1 pgreenhouse1=507/33 logppres1=math.log(ppres1) logpgreenhouse1=math.log(pgreenhouse1) pregrek1=logpgreenhouse1/logppres1 tempk1=737/288 ## temperature #tempk1=288/737 #pmars1=0.0065 #pmars1=92 greenhouse_mars1=33*tempk1*math.exp(pregrek1*math.log(pressu1))

print(pressu1,greenhouse_mars1 )

return(greenhouse_mars1)

def temperature_latitude_in_mars(la1):

TLAT=202.888-0.781801*la1-22.3673*la1*la1-3.16594*la1*la1*la1 return(TLAT)

def calc_tna(Sin): ## NOK

stefan1=5.670374419e-8 ## !!! ne1, ae1 NOK not ok!!!! em1=0.98 ne1=1-0.754 ae1=0.5 a1=0.932*math.pow(ne1, 0.25) a2=math.pow((1-ne1),0.25) a3=2/5*math.pow(((Sin*(1-ae1))/(em1*stefan1)), 0.25) Tna=a3*(a1+a2) return(Tna)

1. 1. attempt to estimete Earhl like atmosph greenhouse temperature

def greenhouse_tna(Sin, pressure): #Sin=1 #pressure= pressure2=pressure*170 ## fre pressure from hat! Tna = 32.44*math.pow((Sin*1367),1/4) Tna2=calc_tna(Sin) #Tna=255 #pcoeff_earth1=1.459 ## for earth 1e5 Pa #pcoeff_mars1=1.194 ± 0.004 # P P # a1=0.174205*math.pow(pressure2, 0.150263) a2=1.83121e-5*math.pow(pressure2, 1.04193) tstna1=math.exp((a1+a2)) Ts=Tna*tstna1 print(" Tna Tna2, Ts", Tna,Tna2, Ts) return(Ts, Tna)

def estimate_greenhouse_temperature_excess(pressure, proportion): ## basis: Earth basic greenhouse ## partially from hat! ## earth co2 2x, T +5.5K ?? ## earth co2 400e-6, 14.7 300e-6 13.7 ## propok=math.log(400/300) ## co2 2x --> T 1.5-4,5 C ## co2 2x, T++ 3 C ## simulated 1d simu ' #PROPTD=math.log(proportion/355e-6)*13.18 ## 329 k moist greenhouse ??? LOG210=3.32192809489 PROPTD=math.log(proportion/355e-6)*LOG210*3 ## 1.5 ...4.5 PROPTD=math.log(proportion/355e-6)*5 ## mean ## 10x 3, 5, 7 id p=1bar ##PROPTD=math.log(proportion/355e-6)*20 ## maybe Earth

# TGH=33+math.log10(pressure)*35 ##*10, 20, 35 K+ if pco2=355 ppm 130 abs max TGH=33+math.log10(pressure)*35+PROPTD ##*10, 20, 35 K+ if pco2=355 ppm 130 abs maxfor earth p 10x t 100 K++ is real! ## for earth p 2x --- t++ 95 K++ moist greenhouse ## earth 1.5 bar, 400 ppm -- T ++ 2K!! 1.7 bar ---> runaway! T++ 100+ K TGHIN=33+math.log10(pressure)*10 TGHMAX=33+math.log10(pressure)*130 ## earth now pressure 0.01 ,0.1, 1, 10, 100 bar ## T 275, 310, 340,410, 540 K ## Runaway limit 0.1 bar atm trunaway 325, 1 bar 360, 10 bar 420, 100 bar 590 K ## runaway net insolation+ground heat ca 285, now it is in earth ca 240 ## earth runaway limit is 320 w/m2 top of atmosphere flux, that is 325 K ## ## p2x 10...20 K+ ## hignest stable influx is 385 w m2 to earthm, S0=1540, that is 1.131 -- >0.94 au

return(TGH)

def calcu_temperature(solarconstant, albedo, emission_ir):

   # tee1 = 394 * math.pow( (1-A),0.25)*math.pow(r1, -0.5)
   stefan = 5.670374419e-8
   # tee1=pow  ( ( ((1-albedo)*fsun) / (4*stefan)) *emission_ir  , 0.25 )
   tee1 = np.power(((solarconstant*(1-albedo)) / (4*stefan*emission_ir)), 1/4)
   return(tee1)

def radiogenic_volcanism_earthlike(mass, radius, agegyr, qratio): ## mercury original heat #https://arxiv.org/ftp/arxiv/papers/1808/1808.01333.pdf

print("Agegyr ", agegyr)

Q0=3.6e-11## W/kg EH condrite model lam=1.5e-17 ## #t=year1*1e9*agegyr #QTOT0=mass*mearth*Q0 ## same than RH #print("QTOT0 TW", QTOT0/1e12) #QTOT=QTOT0*math.exp(-lam*t) area=math.pow((radius*rearth),2)*math.pi*4 ## "empirical" from Earth rh_rel=math.pow(mass, 1.2089) rh_0=197.8*rh_rel*qratio print("radiogenic heat rel ",rh_rel) print(rh_0) tee=year1*agegyr*1e9 rh=rh_0*math.exp(-1*lam*tee) th=rh*2 heatflux=(th*1e12)/area print(" radiogenic heatflux ", heatflux) volcanism=1 ## large volcanism surely off tectonism=1 if(heatflux<0.085): volcanism=0 #https://royalsocietypublishing.org/doi/10.1098/rsta.2017.0416 if(heatflux<0.040): tectonism=0 ## stagnant lid if(heatflux>0.070): tectonism=2 ## moving lid if(heatflux<0.0068): print("Tectonism surely ceased, if not tidal") if(heatflux>0.187): print("Hot stagnant lid possible") if(heatflux>0.482): print("Hot stagnant lid more possible") print("tectonism by heatflux ", tectonism) return(volcanism)

def volcanic_internal_energy(mass): ## radiogenic heating=math.pow(mass,1) ## magma spreading: viscosity, gravity ## high silica spreads wide ? ## mantle size! radius=math.pow(mass, 2.7) #roo=mass/math.pow(radius,2.7) ## plate speed like #vei=math.pow((mass/radius)) ## surface area/vol area 1e-5 htotal 1e10 1e-6 htotal 1e12 area=math.pow(radius,2) volume=math.pow(radius,3) av=area/volume vei_min=1/(av*av) return(vei_min)

def volcanic_ability_index(radius, mass): vai=0 roo=mass/math.pow(radius,2.7) ## plate speed like

##vai=mass/math.pow(radius,2) vai=math.pow(mass, 1.2) return(vai)

def atmosphere_mass_fraction_earthlike(massstar,lstar, massplanet, rau2, r2, agegyr): ## NOK #https://academic.oup.com/mnras/article/504/2/2034/6195519#update1617641408601 rpl=r2+(r2/100) rc=r2 age=agegyr*1e9 ro=1 ## near solar mass 1 ... 1.5

C1=0.0464*math.pow( (lstar/(rau2*rau2)), -0.09)*math.pow(ro, 0.04)*math.pow(massplanet, 0.07 )*math.pow((age/(100*1e6)), 0.03) C2=(1.0296-0.1676)*math.pow(massplanet, 0.095 )*math.pow(ro, 0.01)*math.pow( (lstar/(rau2*rau2)), 0.01) C3=(1.0053-0.1753)*math.pow(massplanet, -0.04 )*math.pow(ro, 0.02)*math.pow( (lstar/(rau2*rau2)), -0.02)*math.pow((age/(100*1e6)), -0.017)

atm_fraction=C1*math.pow((rpl-rc), C2)*math.tanh(C3*(rpl-rc)) print("fraction ",atm_fraction) return(atm_fraction)

def scaleheight(mass, radius, temp, molmass): gee=9.81*(mass/math.pow(radius, 2)) scaleheight=(kb*temp)/(molmass*amu*gee) return(scaleheight);

def radius_planet_50pros_water(mass): # https://people.earth.yale.edu/sites/default/files/files/duffy_madhu_kkml_superEarths_2015.pdf # 50 & rock, 50 %h2O Sotin et al, 2007 radius_50_pros_h2o=1.262*math.pow(mass, 0.275) return(radius)

def color_temperature_to_rgb(kelvin): #https://gist.github.com/paulkaplan/5184275 temp = kelvin/100 red=1 green=1 blue=1 if(temp<=66): red=255 green=temp green=99.4708025861*math.log(green)-161.1195681661

if(temp<=19): blue=0 else: blue=temp-10 blue=138.5177312231*math.log(blue)-305.0447927307 else: red= temp-60 red=329.698727446*math.pow(red,-0.1332047592) green=temp-60 green=288.1221695283*math.pow(green,-0.0755148492) blue=255

r=clamp(red, 0, 255) g=clamp(green, 0, 255) b=clamp(blue, 0, 255) print("Color of star rgb ", r,g,b) return(r,g,b)

def tidal_force_star_planet(mstar1, mass1, distance, radius): tidal_force=(2*G*mstar1*msun*mplanet1*mearth1*distance*au*radius*rearth)/(math.pow((distance*au), 3)) return(tidal_force)

1. def planet_rotation_ratex(mass, agegyr):
2. ## earth original rotation rate 5h, with moon braking
3. ## mass
4. rotrate=math.exp(-agegyr/0.5)*5
5. return(5)

def estimate_planet_rotation_from_mass_and_age(mass1, age1): # ## earth original rotation rate 5h, with moon braking

   # planet rotation by mass, age: no central star tidal braking or collisions etc assumed
   # earth rotation slowing rate
   ## 620 ma ago, 4.0 gyr 21.9+-0.4 h
   ## if earth has no moon, its rotation rate now can meybe be 6-10 h, 6-12 h
   ## tide sun=0,46*tide_moon
   ## 50 gyr earth rotation=31 d ?
   ## 1, 5 gyr, earth rotation rate 0.4 -0.5  x current
   slowratenow = 0.0017/100
   ga = 1000*1000*1000
   # earth rotation rate at beginning
   ##rota1 = 2.448
   rota1=5
   # mass1=15.536
   if (mass1 < 10):
       logmass1 = math.log(mass1)+2
   else:
       logmass1 = math.log(mass1)

   #gerr1 = 2/logmass1
   gerr1 = 1/logmass1*1.748
       
   # kerr1=0.689
   print(age1)

   rotation1 = rota1+slowratenow*gerr1*age1*ga/3600
   rotation2=math.exp(age1/2.9326)*rota1 ## w /moon
   rotation_without_moon1=math.exp(age1/25)*rota1
   rotation_without_moon2=math.exp(age1/5.25)*rota1
   #rotation2=math.exp(age1/2.71)*rota1
   # planet mass in earth masses
   print("rotation2 ", rotation2," h ", rotation2/24, " d")
   print("Rotation, if no moon maybe hours ", round(rotation_without_moon1,4),round(rotation_without_moon2,4))
   return(rotation2)

def calc_tzams_myr(mstar1): tzams1=80.5/math.pow(mstar1, 2.22) return(tzams1)

def calc_lstar_from_mass(mstar1):

lstar1=math.pow(mstar1, 3.5)

if(mstar1>55): lstar1=32000*mstar1 if(mstar1<56): lstar1=1.4*math.pow(mstar1, 3.5) if(mstar1<2): lstar1=1.4*math.pow(mstar1, 4) if(mstar1<0.43): lstar1=0.23*math.pow(mstar1, 2.3)

if(mstar1<0.46): if(mstar1>0.178): lstar1=math.pow(10,(2.028*math.log10(mstar1)-0.976) ) if(mstar1<0.72): if(mstar1>0.45): lstar1=math.pow(10,(4.572*math.log10(mstar1)-0.102) ) if(mstar1<1.05): if(mstar1>0.71): lstar1=math.pow(10,(5.743*math.log10(mstar1)-0.007) ) if(mstar1<2.5): if(mstar1>1.04): lstar1=math.pow(10,(4.329*math.log10(mstar1)+0.010) ) if(mstar1<7.1): if(mstar1>2.4): lstar1=math.pow(10,(3.967*math.log10(mstar1)+0.093) ) if(mstar1<32): if(mstar1>6.9): lstar1=math.pow(10,(2.865*math.log10(mstar1)+1.105) )

return(lstar1)

def read_command_line(): ## command line seed1=12345 seed1=None name1="Planets" outfilename1="planets_0000.csv" type1=1 ## planet system type 1 rocky and icy/gas 2: same wet/icy/gassy mstar1=1 ## 1 sun surf1=1.0 ## 1 snormal metallicity1=1 ## 1 sun alpha1=3e-3 ## tex 1e-2, 1e-3, 1e-4 viscosity alpha for giant planets only! gastau1=15 ## ca 3 normal 1 small plannets gasdep1=7 gastime1=5 ## ca 3 normal 1 small planets migratek1=1 ## default migration_: no migration

first_name = None last_name = None argv = sys.argv[1:] try: opts, args = getopt.getopt(argv, "m:s:") except: print("Error")

for opt, arg in opts: if opt in ['-s']: seed1 = int(arg) elif opt in ['-m']: mstar1 = float(arg) elif opt in ['-M']: metallicity1 = float(arg) elif opt in ['-D']: surf1 = float(arg) elif opt in ['-t']: type1 = int(arg) elif opt in ['-time']: gastime1 = float(arg) elif opt in ['-taugas']: gastau1= float(arg) elif opt in ['-taudep']: deptau1= float(arg) elif opt in ['-alpha']: alpha1= float(arg) elif opt in ['-migratek']: migratek1= float(arg) elif opt in ['-drawskala']: drawskala1= float(arg) elif opt in ['-N']: name1 = arg elif opt in ['-o']: outfilename1, name1 = arg

#print(name1, outfilename1, type1, seed1, mstar1, metallicity1, surf1,gastime1, gastau1, alpha1, migratek1)

return(name1, outfilename1, type1, seed1, mstar1, metallicity1, surf1,gastime1, gastau1,deptau1, alpha1, migratek1, drawskala1)

def gravities_earths(m1, r1): return(m1/(r1*r1))

def vescs_earths(m1, r1): return(np.sqrt(m1/r1))

def rotation_periods_hours(mass, radius): ## attempt to empirical formula rotdays=1*np.power(mass,-1/8.3) return(rotdays*24)

def rotation_periods_hours_teor(mass, radius): # https://arxiv.org/pdf/2301.13836.pdf rotdays=0.632*np.power(mass, 1/2)*radius # rotdays= 0.632*math.pow(mass, 1/6)*math.pow(rooearths, 1/3) return(rotdays*24)

def calculate_tidal_heat_1(mcentral1, asat1, msat1,rsat1, eccsat1): meanmot1=math.sqrt((mcentral1+msat1)/(asat1*asat1*asat1)) etidal_rel1=(mcentral1*mcentral1)*math.pow(rsat1, 5)*meanmot1*eccsat1*eccsat1/math.pow(asat1,6) etidal_rel2=etidal_rel1/6.457868774140072e-16 ## io unit, stony volcanism ## 1 io unit stony volcanism ## 0.2 io unit europa: past stony, current ice volcanism ## 0.00048 enceladus ice volcanism and craters # 0.00061 rhea some cracks ## 2.65 e-5 titania return(etidal_rel2)

def current_milli_time():

   return round(time.time() * 1000)

def radius_s1(ms1, rockp1, icep1, gasp1): #rr1=math.pow((ms1*rockp1/100), 0.27) ## colod hygrogen earth radius 4 if(ms1>(80*320)): ms1=80*320 rr1=math.pow(abs(ms1), 0.27) ## o.4 to 1, not 0.4 per cent if(icep1>0.4): rr1=rr1*1.5 ## coarse estim

if(ms1>5):rr1=math.pow((ms1), 0.55)*0.75 ## density below 3.3 exponent 0.63 ? if(ms1>50):rr1=math.pow((ms1), 0.02)*7 ## 8 ? if(ms1>1500):rr1=math.pow((ms1/300), -1/8)*11.2 #density1=rockp1/100+(icep1/100)*0.2+(gasp1/100)*0.05 #volume1=ms1/density1 #print(density1) #rr1=rr1*math.pow(volume1, 0.2) return(rr1)

def estimate_planet_density_from_distance_from_sun(rau1):

   # planet density, f distance of sun
   roo1 = 4100*math.pow(rau1, -0.21)  # R2=0.88
   return(roo1)

def can_hold_some_time_gas_min_molmass(mass, radius, temp): mmolmin=54*temp*(kb*radius*rearth)/(mass*mearth*G) mmolmin=mmolmin/amu return(mmolmin)

def can_hold_long_time_gas_min_molmass(mass, radius, temp): mmolmin=400*temp*(kb*radius*rearth)/(mass*mearth*G) mmolmin=mmolmin/amu return(mmolmin)

def calculate_radius_from_mass_stony_planet(mass): radius=1*math.pow(mass, 0.2739-(0.008*math.log(mass)) ) return(radius)

def uncompressed_planet_material_density(lumasol, distanceau): arel=math.sqrt(lumasol)*math.sqrt(distanceau) ## estimated uncompressed, unporous planetesimal density roo=4100*math.pow(arel, -0.21)/rooearth return(roo)

def is_tidal_locked(agegyr,massstar,massplanet, rau2, rotperiod): #https://worldbuilding-workshop.fandom.com/wiki/Worldbuilding_Guide_Part_7:_Atmospheres_%26_Oceans?mobile-app=true&theme=dark' # https://www.as.utexas.edu/astronomy/education/spring02/scalo/kasting.pdf p1=np.sqrt(((rau2*rau2*rau2)/massstar)) #Q1=100 ## earth today 13, can be 100 #Q1=13 Q1=10

rt=0.027*np.power( ((p1*agegyr*1e9)/Q1),1/6)*math.pow(massstar, 1/3) print("rt ", rt) #ee=(massstar*agegyr)/massplanet #locked=0

#if(rau2<rt): locked=1

locked=np.copy(rt)

len1=len(locked)

for n in range(0,len1): if (locked[n]>rau2[n]): locked[n]=1 else: locked[n]=0

return(locked)

def minimum_molecule_mass_to_hold_in_atmosphere(temperature, mass, radius): molekm=(150*kb*temperature*radius*rearth)/(G*mass*mearth) ## atmosphere minimum pressure, with oxygen mask 0.145 atm, in 13700 meters return(molekm/amu)

def calculate_planet_density_radius_from_mass_amounts_1(fe0, stone0,ice0,water0, gas0): iron_density1=7.79 stone_density1=4.54 ice_density1=1.5 gas_density1=0.1 water_density1=1 mass1=fe0+stone0+ice0+gas0 fe1=fe0/mass1 stone1=stone0/mass1 ice1=ice0/mass1 gas1=gas0/mass1 water1=water0/mass1

density1=fe1*iron_density1+stone1*stone_density1+ice1*ice_density1+water_density1*water1+gas1*gas_density1 radius1=math.pow((mass1/(density1/5.51)), 1/3)*rearth

return(mass1*mearth, density1*1000, radius1)

def magnetic_potential_estimation(mass, radius, rotperiod):

   density=mass/(radius*radius*radius)
   omega=24/rotperiod
   #alfa=0.15 ## dynamo 0.05 multipolar other 1
   alfa=1
   magnetic=np.power(density, 1/2)*np.power(radius,3)*alfa *np.power(mass, 1/8)*omega## math pow mass, 1/6 is  omega, rotaion angle freq
   return(magnetic)

def atmosphere_pressure_estimation(mass, radius, temp): gravity=mass/np.power(radius,2) pressure = gravity*mass*radius*radius #*np.sqrt(temp/288) return(pressure)

def calc_esi_term(o,r,we): w=we/4 a=abs((o-r)/(o+r)) valu=np.power((1-a),w) return(valu)

def calc_esis(mass, radius, temp):

w_radius=0.57 w_density=1.07 w_vesc=0.7 w_temp=5.58 ref_temp=288

vesc=np.sqrt(mass/radius) density=mass/(radius*radius*radius)

esi_i=calc_esi_term(radius,1,w_radius)*calc_esi_term(density,1,w_density) esi_s=calc_esi_term(vesc,1,w_vesc)*calc_esi_term(temp,ref_temp,w_temp)

esi=esi_i*esi_s

#return(esi, esi_i,esi_s) return(esi)

def planet_temperatures_1(starmass1, starlum1, aus1):

   ## temp, Kelvin
   aus2=aus1*np.power(starlum1, -1/2)*math.pow(starmass1, -1/7) ## check mass effect!
   temps1=np.power(aus2, -1/2)*288
   return(temps1)

def increase_mass_from_dustgas_1(mstar1,snowline1, incdustcoeff1, incgascoeff1, bhill1, time1, timestep, taugas, tausolids, aas1, eccs1, ms1, mrocks1, mices1, mgases1, diska0, dust1, gas1, ices1, rocks1):

aas2, eccs2, ms2, mrocks2, mices2, mgases2, dust2, gas2, ices2, rocks2=accrete_mass_from_dustgas_1(mstar1,snowline1, incdustcoeff1, incgascoeff1, bhill1, time1, timestep, taugas, tausolids, aas1, eccs1, ms1, mrocks1, mices1, mgases1, diska0, dust1, gas1, ices1, rocks1)

return(aas2, eccs2, ms2, mrocks2, mices2, mgases2, dust2, gas2, ices2, rocks2)

def combine_by_relative_distance_base_1(aas1, eccs1, ms1,mrocks1, mices1, mgases1, idx1, idx2, dist1): #print(len(aas1), idx1, idx2) #print("LL") #print (len(aas1), len(eccs1), len(ms1), len(mrocks1), len(mices1), len(mgases1)) df1=pd.DataFrame(data=np.dstack([aas1, eccs1, ms1, mrocks1, mices1, mgases1])[0] ) #print("LLL") #print(df1)

df2=df1.sort_values(by=0)

#print(df1) #quit(-1)

if(idx1==idx2): #print("samedex") return (aas1, eccs1, ms1, mrocks1, mices1, mgases1)

if(idx1>idx2): idx0=idx1 idx1=idx2 idx2=idx0

a1=df2[0][idx1] a2=df2[0][idx2] e1=df2[1][idx1] e2=df2[1][idx2] m1=df2[2][idx1] m2=df2[2][idx2] mrock1=df2[3][idx1] mrock2=df2[3][idx2] mice1=df2[4][idx1] mice2=df2[4][idx2] mgas1=df2[5][idx1] mgas2=df2[5][idx2]

if(a1>a2): idx0=idx1 idx1=idx2 idx2=idx0 a1=df2[0][idx1] a2=df2[0][idx2] e1=df2[1][idx1] e2=df2[1][idx2] m1=df2[2][idx1] m2=df2[2][idx2] mrock1=df2[3][idx1] mrock2=df2[3][idx2] mice1=df2[4][idx1] mice2=df2[4][idx2] mgas1=df2[5][idx1] mgas2=df2[5][idx2]

#print("a1", a1) #print("a2", a2)

diff1=abs(a2-a1)

am1=(a1+a2)/3 dist2=am1*dist1

if(diff1<dist2): diffa1=a2-a1 a3=(a1+a2)/2 m3=m1+m2 mrock3=mrock1+mrock2 mice3=mice1+mice2 mgas3=mgas1+mgas2

e3=(e1+e2)/2 ms1=m1/m3 ms2=m2/m3

a3=a1+ms1*diffa1 e3=e1*ms1+e2*ms2

df2[0][idx1]=a3 df2[1][idx1]=e3 df2[2][idx1]=m3 df2[3][idx1]=mrock3 df2[4][idx1]=mice3 df2[5][idx1]=mgas3 df3=df2.drop(df2.index[idx2]) else: df3=df2

aas1=np.array(df3[0]) eccs1=np.array(df3[1]) ms1=np.array(df3[2]) mrocks1=np.array(df3[3]) mices1=np.array(df3[4]) mgases1=np.array(df3[5]) #print("--->") #print(aas1) #print(eccs1) #print(ms1)

#quit(-1) return (aas1, eccs1, ms1, mrocks1, mices1, mgases1)

def combine_by_relative_distance_1(aas1, eccs1, ms1, mrocks1, mices1, mgases1,dist1): len1=len(aas1) nummula=int(len1*len1/2) for n in range (0, nummula): len1=len(aas1) idx1=np.random.randint((len1-1)) idx2=np.random.randint((len1-1))

aas2, eccs2, ms2, mrocks2, mices2, mgases2=combine_by_relative_distance_base_1(aas1, eccs1, ms1, mrocks1, mices1, mgases1,idx1, idx2, dist1) aas1=np.copy(aas2) eccs1=np.copy(eccs2) ms1=np.copy(ms2) mrocks1=np.copy(mrocks2) mices1=np.copy(mices2) mgases1=np.copy(mgases2) #print (len(aas2), len(eccs2), len(ms2), len(mrocks2), len(mices2), len(mgases2))

#print("MAMMAGAMMA") #print (len(aas1), len(eccs1), len(ms1), len(mrocks1), len(mices1), len(mgases1))

return (aas1, eccs1, ms1, mrocks1, mices1, mgases1)

def find_nearest_resonances_by_external_1(aas1,ms1, aext1): ## NOK ## attempt to adjust distances of born planets by resonances planetab1=np.dstack((aas1, ms1)) planetab2=planetab1[0]

#print(planetab2) #planetab2=planetab2[planetab2[:,1].argsort()] #planetab2=np.sort(planetab2, order=['f1'], axis=1) #planetab3=planetab2[np.lexsort(planetab2)] #print(planetab3)

#df1 = pd.DataFrame(planetab2, index=planetab2[:,0]) df1 = pd.DataFrame(planetab2) df1.colums=["aas", "ms"] df1=df1.sort_values(0) #print(df1) #quit(-1)

aas1=np.array(df1[0]) ms1=np.array(df1[1])

#print(aas1) #print(ms1) #quit(-1)

#aas1=planetab3[:,0] #ms1=planetab3[:,1]

aas2=np.copy(aas1)*0 len1=len(aas1) firstreso1=1/48 abase1=aext1*math.pow(firstreso1, 2/3) #print(abase1) reso1=firstreso1 #for n in range(0,len1 ): # a1=aext1*math.pow(reso1, 2/3) # aas2[n]=a1 # reso1=reso1*2

## second try

for n in range(0,len1-1 ): a1=aas1[n] a2=aas1[n+1] m1=ms1[n] m2=ms1[n+1] pm1=m1 pm2=m2 if(pm2<pm1): pm3=pm1 pm1=pm2 pm2=pm3 mrel1=pm2/pm1 #print("*",mrel1,pm1,pm2) ta1=aext1*math.pow(reso1, 2/3) resok1=1.8 if(mrel1<1.3): resok1=1.3 if(mrel1>2): resok1=2.0 if(mrel1>4): resok1=2.1 if(mrel1>8): resok1=2.2 if(mrel1>10000): mrel1=2 reso1=reso1*resok1 ta2=aext1*math.pow(reso1, 2/3)

aas2[n]=ta1 aas2[n+1]=ta2

#print(aas2) #quit(-1) return(aas2)

def consentrate_dust_to_narrow_rings_for_making_planets(adisk1, dust1, rocks1, ices1):

num1=len(dust1) deltaecc1=0.2 amax1=max(adisk1) amin1=min(adisk1) adelta1=adisk1[1]-adisk1[0] awidth1=amax1-amin1 amean1=(amax1+amin1)/2

for i in range(0,1): #print("I ", i) dust2=np.copy(dust1) rocks2=np.copy(rocks1) ices2=np.copy(ices2)

ranke1=int((amax1-amin1)/(deltaecc1*amean1)) ranke1=ranke1*2 #ranke1=1 for m in range(0, ranke1): p1=np.random.randint(num1) #print(p1) #quit(-1) #p1=int(np.random.random()*(num1/2))+int(num1/2) #p1=int(np.random.normal(120,25)) if(p1<0): p1=0 if(p1>num1): p1=num1-1 #delta1=int(p1*deltaecc1)*2 delta1=int((deltaecc1*num1)/2) p2=p1-delta1 p3=p1+delta1 if(p2<0): p2=0 if(p3>num1): p3=num1-1 #print("delta1 ", delta1) #print("p ", p1, p2, p3) mass1=0 mice1=0 mrock1=0 for n in range(p2, p3, 1): dusta1=dust1[n] ice1=ices1[n] rock1=rocks1[n] dustb1=dusta1/2 dustc1=dusta1-dustb1 mass1=mass1+dusta1 mice1=mice1+ice1 mrock1=mrock1+rock1 dustc1=0 dust1[n]=0 dust2[n]=0

dust2[p1]=mass1 rocks2[p1]=mrock1 ices2[p1]=mice1

dust1=np.copy(dust2) ices1=np.copy(ices2) rocks1=np.copy(rocks2)

return(dust2, rocks2, ices2)

def collect_dust_to_planetesimals_lane_by_lane(snowline1, prop1, adisk1, dust1, rock1, ice1):

aas1=[] ms1=[] mrocks1=[] mices1=[]

len1=len(adisk1)

for n in range(0, len1): a1=adisk1[n] m1=dust1[n]*prop1

icep1=0.5 rockp1=0.5

if(a1<snowline1): rockp=1 icep1=0

mice1=m1*icep1 mrock1=m1*rockp1

dust1[n]=0 if(m1>0): aas1.append(a1) ms1.append(m1) mrocks1.append(mrock1) mices1.append(mice1)

aas1=np.array(aas1) ms1=np.array(ms1) mrocks1=np.array(mrocks1) mices1=np.array(mices1)

return(dust1, aas1, ms1, mrocks1, mices1)

def explode_planetesimal_to_dust1(adisk1, snowline1, dust1, drock1, dice1, a1, m1,mrock1, mice1, dustecc1): len1=len(dust1) ai1=np.searchsorted(adisk1, a1) #print(ai1) #dust1[ai1]=m1 dustk1=(m1*2)/len1 dustadd1=norm.pdf(adisk1, loc=a1, scale =dustecc1)*dustk1 darock1=np.copy(dustadd1) daice1=np.copy(dustadd1) mtotal1=mrock1+mice1 prock1=mrock1/mtotal1 pice1=dice1/mtotal1 darock1=darock1*prock1 daice1=daice1*pice1 icek1=np.copy(adisk1)

icek1=np.where(icek1>snowline1,0.0,1 ) iceadd1=daice1*icek1 rockadd1=deice1*icek1 ## assumes ices in snowline to gas !!! #dust1=dust1+dustadd1

dust1=dust1+dustadd1 drock1=drock1+rockadd1 dice1=dice1+iceadd1

sum1=np.sum(dust1) print(sum1) return(adisk1, dust1, drock1, dice1)

def collect_nearby_planetesimals_from_ecc_and_hill_radius(random1, mstar1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, idx1, bee1, adisk1, dust1, ices1, rocks1): ## NOTE *** this is basic combination of planetesimal routine mstar1_me=333000*mstar1 aas2=[] ecc2=[] ms2=[] mrocks2=[] mices2=[] mgases2=[] maksi1=len(aas1)-1

select1=[]

if (maksi1==0): return(aas1, eccs1, ms1, mrocks1, mices1, mgases1, dust1, ices1, rocks1)

for n in range(0, maksi1): #if(idx1!=n): a1=aas1[idx1] e1=eccs1[idx1] m1=ms1[idx1] mrock1=mrocks1[idx1] mice1=mices1[idx1] mgas1=mgases1[idx1]

a2=aas1[n] e2=eccs1[n] m2=ms1[n] mrock2=mrocks1[n] mice2=mices1[n] mgas2=mgases1[n]

if(a1>a2): a3=a1 e3=e1 m3=m1 a1=a2 e1=e2 m1=m2 a2=a3 e2=e3 m2=m3

deltaa1=e1*a1 deltaa2=e2*a2

rhill1=a1*math.pow((m1/(3*(m1+mstar1_me))),1/3 ) addhill1=rhill1*bee1 rhill2=a2*math.pow((m2/(3*(m2+mstar1_me))),1/3 ) addhill2=rhill2*bee1

addhills12=addhill1+addhill2

aapo1=a1+deltaa1+addhills12 aperi2=a2-deltaa1-addhills12

#aapo1=a1+deltaa1 #aperi2=a2-deltaa1

randi1=np.random.randint(10000) ## not all tries to combine ... ## randi 10 #if(a2!=a1): #if(n!=idx1): if(random1==0): ## append always if(aperi2<aapo1): select1.append(n) else: ## append by randoms: eq collision and combining a ## more rare than always if(randi1<5): ## eq 10 if(aperi2<aapo1): select1.append(n)

#print("select :", select1)

#quit(-1)

selen1=len(select1) if (selen1==0): #print("selen 0") return(aas1, eccs1, ms1, mrocks1, mices1, mgases1, dust1, ices1, rocks1)

aa1=0 mm1=0 ee1=0

## "central object" a0=aas1[idx1] e0=eccs1[idx1] m0=ms1[idx1] mrock0=mrocks1[idx1] mice0=mices1[idx1] mgas0=mgases1[idx1]

at1=0 et1=0 mt1=0 mrockt1=0 micet1=0 mgast1=0

tolen1=0 #print("summing ", selen1)

for n in range(0, selen1): #print(idx1, select1[n]) sedex1=select1[n] if(idx1!=n): a1=aas1[sedex1] e1=eccs1[sedex1] m1=ms1[sedex1] mrock1=mrocks1[sedex1] mice1=mices1[sedex1] mgas1=mgases1[sedex1] at1=at1+a1 et1=et1+e1 mt1=mt1+m1 mrockt1=mrockt1+mrock1 micet1=micet1+mice1 mgast1=mgast1+mgas1 tolen1=tolen1+1

#print("sum mt1", mt1) at1=at1/selen1 et1=et1/selen1 ## WARNING no orbit circulation! mt1=mt1 mrockt1=mrockt1 micet1=micet1 mgast1=mgast1

tta1=a0 tte1=e0 ttm1=m0 ttmrock1=mrock0 ttmice1=mice0 ttmgas1=mgas0

#quit(-1)

for n in range(0, selen1): #print(select1[n]) #print(sedex1, n) if(idx1!=n): sedex1=select1[n] wa1=aas1[sedex1] we1=eccs1[sedex1] wm1=ms1[sedex1] wmrock1=mrocks1[sedex1] wmice1=mices1[sedex1] wmgas1=mgases1[sedex1] da1=a0-wa1 de1=e0-we1 mk1=wm1/mt1 deltawa1=wa1-a0 absdeltawa1=abs(deltawa1) sign1=0

if(deltawa1!=0): sign1=(wa1-a0)/(abs(wa1-a0))

masak1=m0/(wm1+m0)

tta1=tta1+mk1*da1*sign1*masak1 tte1=tte1+mk1*de1*sign1*masak1 ttm1=ttm1+wm1 ttmrock1=ttmrock1+wmrock1 ttmice1=ttmice1+wmice1 ttmgas1=ttmgas1+wmgas1

at1=at1+a1 et1=et1+e1 mt1=mt1+m1 mrockt1=mrockt1+mrock1 micet1=micet1+mice1 mgast1=mgast1+mgas1

#quit(-1) for n in range(0, selen1): #print(select1[n]) sedex1=select1[n] aas1[sedex1]=0 eccs1[sedex1]=0 ms1[sedex1]=0 mgases1[sedex1]=0 mrocks1[sedex1]=0 mices1[sedex1]=0

#tta1=at1/selen1 #tte1=et1/selen1 #ttm1=mt1 #print(tta1, ttm1, mt1)

#at1=at1/selen1 #et1=et1/selen1 ## WARNING no orbit circulation! #mt1=mt1

#print("*") #print(at1) #print(et1) #print(mt1)

#print(tta1) #print(tte1) #print(ttm1)

## "central object" aas1[idx1]=tta1 eccs1[idx1]=tte1 ms1[idx1]=ttm1 mrocks1[idx1]=ttmrock1 mices1[idx1]=ttmice1 mgases1[idx1]=ttmgas1

#quit(-1)

for n in range(0, (maksi1+1)): #print(n) a1=aas1[n] e1=eccs1[n] m1=ms1[n] mrock1=mrocks1[n] mice1=mices1[n] mgas1=mgases1[n]

if(a1==0): a1=00 else: aas2.append(a1) ecc2.append(e1) ms2.append(m1) mrocks2.append(mrock1) mices2.append(mice1) mgases2.append(mgas1)

aas3=np.array(aas2) ecc3=np.array(ecc2) ms3=np.array(ms2) mrocks3=np.array(mrocks2) mices3=np.array(mices2) mgases3=np.array(mgases2)

#print(aas3) #print(ms3) #print(mrocks3) #print(mices3) #print(mgases3) #quit(-1)

return(aas3, ecc3, ms3, mrocks3, mices3, mgases3, dust1, ices1, rocks1)

def combine_small_planetesimals_to_bigger_objects(aas1, eccs1, ms1,mrocks1, mices1, mgases1, mlimit1, adisk1, dust1, ices1, rocks1): maksi1=len(aas1)

select1=[]

for n in range(0, maksi1): a1=aas1[n] e1=eccs1[n] m1=ms1[n]

if(m1>mlimit1): select1.append(n)

#print(select1) lenselect1=len(select1)

pas1=[] peccs1=[] pms1=[] pmrocks1=[] pmices1=[] pmgases1=[]

pas1.append(0) peccs1.append(0) pms1.append(1e-6) ## inject small planetesimal to pmrocks1.append(1e-6) ## inject small planetesimal to pmices1.append(0.0) ## inject small planetesimal to pmgases1.append(0.0) ## inject small planetesimal to # make algo to function

for m in range(0, lenselect1): idx1=select1[m] a1=aas1[idx1] e1=eccs1[idx1] m1=ms1[idx1] mice1=mices1[idx1] mrock1=mrocks1[idx1] mgas1=mgases1[idx1] #print(a1, e1, m1) pas1.append(a1) peccs1.append(e1) pms1.append(m1) pmrocks1.append(mrock1) pmices1.append(mice1) pmgases1.append(mgas1)

a2=max(np.array(aas1))*1.6 ## assume small planetoid to outsize study area ## to make algo ok! pas1.append(a2) peccs1.append(0.0) pms1.append(1e-5) ## inject small planetesimal to pmrocks1.append(1e-5) pmices1.append(0) pmgases1.append(0)

pas1=np.array(pas1)

peccs1=np.array(peccs1) pms1=np.array(pms1) pmrocks1=np.array(pmrocks1) pmices1=np.array(pmices1) pmgases1=np.array(pmgases1)

plen1=len(pas1)

#print(pas1) #print(pms1)

#quit(-1)

planetab1=np.dstack((pas1, peccs1, pms1, pmrocks1, pmices1, pmgases1)) planetab2=planetab1[0]

planetab2=planetab2[planetab2[:,0].argsort()]

pas1=planetab2[:,0] peccs1=planetab2[:,1] pms1=planetab2[:,2] pmrocks1=planetab2[:,3] pmices1=planetab2[:,4] pmgas1=planetab2[:,5]

#print("¤¤¤") #print(pas1) #print(peccs1) #print(pms1) #quit(-1)

planets1=len(pas1)-1 ## jn warn!

pas2=[] peccs2=[] pms2=[] pmrocks2=[] pmices2=[] pmgases2=[]

#print("**") for m in range(1, (planets1-1)): a1=pas1[m-1] e1=peccs1[m-1] m1=ms1[m-1] mrock1=mrocks1[m-1] mice1=mices1[m-1] mgas1=mgases1[m-1]

a2=pas1[m] e2=peccs1[m] m2=pms1[m] mrock2=mrocks1[m] mice2=mices1[m] mgas2=mgases1[m]

a3=pas1[m+1] e3=peccs1[m+1] #m3=ms1[m+1] #mrock3=mrocks1[m+1] #mice3=mices1[m+1] #mgas3=mgases1[m+1]

ain1=(a1+a2)/2 aout1=(a2+a3)/2 #print(a1, a2, a3) #print("-->", a2, ain1, aout1) muksu1=len(aas1)-1 smalls1=[]

maksi1=len(aas1)-1 for n in range(0, maksi1): a4=aas1[n] e4=eccs1[n] m4=ms1[n] mrock4=mrocks1[n] mice4=mices1[n] mgas4=mgases1[n] #print(a4, ain1, aout1)

ine1=0 disrupt1=0 if(a4>ain1): if(a4<aout1): ine1=1

if(m!=n): if(ine1==1): smalls1.append(n)

smallnum1=len(smalls1) for i in range(0, smallnum1): idx1=smalls1[i] a5=aas1[idx1] e5=eccs1[idx1] m5=ms1[idx1] mrock5=mrocks1[idx1] mice5=mices1[idx1] mgas5=mgases1[idx1]

## do not change bigger planets values other than mass m2=m2+m5 mrock2=mrock2+mrock5 mice2=mice2+mice5 mgas2=mgas2+mgas5

aas1[idx1]==0 eccs1[idx1]=0 ms1[idx1]=0 mrocks1[idx1]=0 mices1[idx1]=0 mgases1[idx1]=0

pas2.append(a2) peccs2.append(e2) pms2.append(m2) pmrocks2.append(mrock2) pmices2.append(mice2) pmgases2.append(mgas2)

pas2=np.array(pas2) peccs2=np.array(peccs2) pms2=np.array(pms2) pmrocks2=np.array(pmrocks2) pmices2=np.array(pmices2) pmgases2=np.array(pmgases2)

#print(pas2) #quit(-1)

return(pas2, peccs2, pms2, pmrocks2, pmices2, pmgases2,dust1, ices1, rocks1)

def plotting4(aas, eccs, ms, x1,x2): len1=len(aas) kkk=200 #print(aas) #print(ms) plt.xlim((x1, x2)) for n in range(0, len1): a1=aas[n] m1=ms[n] facecolor1="green" if(m1<12): facecolor1="lightblue" if(a1<3): facecolor1="red" plt.scatter(a1,0, s=np.power(m1, 1/2)*kkk, alpha=0.5, facecolor=facecolor1, color=facecolor1, linewidth=3)

return(0)

def plotting3(aas, eccs, ms, x1,x2): kkk=200 print(aas) #print(ms) plt.xlim((x1, x2)) #plt.scatter(aas,aas*0+0, s=np.power(ms, 1/2)*kkk, facecolor="none", color="#ff0000", linewidth=5) plt.scatter(aas,aas*0+0, s=np.power(ms, 1/2)*kkk, alpha=0.3, facecolor="green", color="#007f00", linewidth=2) #plt.scatter(aas,aas*0+1, s=ms*kkk)

def plotting2(aas, eccs, ms, x1,x2): kkk=2000 print(aas) #print(ms) plt.xlim((x1, x2)) plt.scatter(aas,aas*0+0, s=np.power(ms, 1/2)*kkk, facecolor="none", color="#ff0000", linewidth=5) #plt.scatter(aas,aas*0+1, s=ms*kkk)

def plotting(aas, eccs, ms, x1,x2): kkk=2000 plt.xlim((x1, x2)) plt.scatter(aas,aas*0+0, s=np.power(ms, 1/2)*kkk, alpha=0.5)

1. 1. seems ok, little inward migration alos becomes ....

def find_nearest_resonances_1(aas1, accuracy1): # aas2=find_nearest_resonances_1(aas1, 1/float(len(aas1)*2)) aas2=np.copy(aas1)*0 len1=len(aas1) aas2[0]=aas1[0] for i in range(1, len1): a1=aas1[i-1] a2=aas1[i] trojan=0 if(a1==a2): trojan=1 if(trojan==0): p1=math.pow(a1, 3/2) p2=math.pow(a2, 3/2) ratio1=p2/p1 matar=4*int(1/(accuracy1)) kliffa=0

for p in range(matar, 1, -1): for q in range(matar, 1, -1): kliffa=0 dreso1=0 samed=0 if(p==q): samed=1 reso=1 if(samed==0): reso1=q/p dreso1=abs(reso1-ratio1) if(dreso1<accuracy1): kliffa=1 break

if(kliffa==1): break

if(kliffa==1): p3=p1*reso1 print("") print(p1,p2,p3) print(p, q, reso1, ratio1, dreso1) aas2[i]=a1*math.pow(reso1, 2/3) if(kliffa==0): aas2[i]=-99

return(aas2)

def migrate_inward2(aas1, ms1, migratek1): migratek2=(ms1*migratek1)/(aas1*aas1) aas2=aas1*migratek2

if(aas2<0.001): aas2=0.001

return(aas2)

def migrate_inward_1(aas1, migratek1): aas2=aas1*migratek1 if(aas2<0.001): aas2=0.001 return(aas2)

def sculpt_disk_by_planets_1(aas1, ms1, diska0, dust1, gas1, ices1, rocks1): ## maybe nok aas2=np.copy(aas1) ms2=np.copy(ms1) dust2=np.copy(dust1) gas2=np.copy(gas1) ices2=np.copy(ices1) rocks2=np.copy(rocks1) mgap1=30.0 len1=len(aas1) #aas1[5]=10 #ms1[5]=100 for n in range(0, len1): a1=aas1[n] m1=ms1[n] ## location in disk array idx1 = (np.abs(diska0 - a1)).argmin() idx2=idx1-1 #print(m1) if(m1>mgap1): print("GAP ", m1) sk1=0.1 gas2[0:idx2]=gas2[0:idx2]*sk1 dust2[0:idx2]=dust2[0:idx2]*sk1 ices2[0:idx2]=ices2[0:idx2]*sk1 rocks2[0:idx2]=rocks2[0:idx2]*sk1

#plt.plot(diska0, gas2) #plt.show() #quit() return(aas2, ms2, dust2, gas2, ices2, rocks2)

def set_planets_resonances_by_biggest_planet_1(aas1, ms1): ## maybe nok aas2=np.copy(aas1)*0 ms2=np.copy(ms1) accuracy1=1/float(len(aas1)*2) indices1 = np.where(ms1 == ms1.max()) print(indices1) idx1=indices1[0]

print("MAX ", idx1, aas1[idx1], ms1[idx1]) len1=len(aas1) aas2[idx1]=aas1[idx1] for i in range(0, len1): a1=aas1[i] a2=aas1[idx1] if(a2<a1): a1=aas1[idx1] a2=aas1[i]

trojan=0 if(a1==a2): trojan=1 if(trojan==0): p1=math.pow(a1, 3/2) p2=math.pow(a2, 3/2)

ratio1=p2/p1 #matar=1*int(1/(accuracy1)) matar=100 kliffa=0

for p in range(2, 48, 1): for q in range(2, 48, 1): kliffa=0 dreso1=0 samed=0 if(p==q): samed=1 reso=1 if(samed==0): reso1=q/p dreso1=abs(reso1-ratio1) if(dreso1<accuracy1): kliffa=1 break

if(kliffa==1): break

if(kliffa==1): p3=p1*reso1 #print("") #print(p1,p2,p3) #print(p, q, reso1, ratio1, dreso1) aas2[i]=a1*math.pow(reso1, 2/3) if(kliffa==0): #aas2[i]=-99 aas2[i]=aas1[i] if(reso==1): aas2[i]=aas1[i]

print(aas1) print(aas2) return(aas2, ms2)

def gasdisk_cavity1_lab1(): ## nok diskbase1=10 beta1=-3/2 #beta1=0 ain1=0.01 aout1=100 perau1=100 slices1=int((aout1-ain1)*100) aas1=np.linspace(ain1, aout1, slices1) sigmagas1=np.power(aas1/diskbase1, beta1) as1=[5] ms1=[300] idx1=0 ias1=int((as1[idx1]-ain1)*perau1) gapk1=1/ms1[0] ## create inner hole ##sigmagas1[0:ias1]=sigmagas1[0:ias1]*gapk1 gapper1=np.copy(aas1)*1 mu, sigma = 1, 0.05 # mean and standard deviation noise1 = np.random.normal(mu, sigma, slices1) #plt.plot(s) gapp1=1-norm.pdf(aas1, loc=5, scale =0.5) cavity1=norm.cdf(aas1, loc=5, scale =0.5) #plt.plot(aas1, cavity1,'r-', lw=5, alpha=0.6, label='norm pdf') plt.show() #quit(-1) sigmagas1=sigmagas1*cavity1*gapp1*noise1 #sigmagas1=sigmagas1*gapp1 plt.plot(aas1, sigmagas1) plt.xlim((0.1, 30)) plt.ylim((0,2)) plt.show() return(0)

def resonance_stripes(): ## simple study of resonnece stripes in dust disk .... ainner1=0.01 aouter1=5 adelta1=0.01 sigma0=1 aas1=np.arange(ainner1, aouter1, adelta1) pinta1=np.copy(aas1)*0 len1=len(aas1) abase1=5.2 for m in range(1,64*2): for n in range(1,64*2): reso1=n/m reso2=m/n areso1=math.pow(reso1, 2/3)*abase1 areso2=math.pow(reso2, 2/3)*abase1 loc1=int((areso1-ainner1)/adelta1) loc2=int((areso2-ainner1)/adelta1) if(loc1<0): loc1=0 if(loc1>(len1-1)): loc1=0 if(loc2<0): loc2=0 if(loc2>(len1-1)): loc2=0 pinta1[loc1]=pinta1[loc1]+1 pinta1[loc2]=pinta1[loc2]+1

pinta1[0]=0

zeros1=np.where(pinta1==0)

print(zeros1)

smalls1=np.where(pinta1<3) print(smalls1) plt.plot(aas1, pinta1) plt.show() return

def generate_inner_planets_1(mstar1, surf1, metallicity1, taugas1, taudep1, lifetime1, migratek1):

incdustcoeff1=1 incgascoeff1=1

maxplanets=8 ## 8 bee1=2.4 bhill1=2.4 #sigmask1=mstar1*surf1*0.003 #sigmask1=mstar1*surf1*70 ## orig sigmask1=mstar1*surf1*60 f=0.013*metallicity1 taugas1=1e5 ## note taugas not effect taudep1=1e5 tausolids1=1e6 ## not used timebegin=0 timemax=int(lifetime1*1e6) timestep=100000

migk1=migratek1 ## migration 1 no migration, 0 falls to star!

mlimit1=0.02 ## min planet mass ainner1=0.4*math.sqrt(mstar1) aouter1=2.2*math.sqrt(mstar1) adelta1=0.01*math.sqrt(mstar1) diskdelta1=0.01*math.sqrt(mstar1) #meanmassk1=2e-7*2 meanmassk1=3e-10 ## 5 e-10 mean mass, stony planetesimals, icy has ca 4x maxmassk1=meanmassk1*1 minmassk1=meanmassk1*1 maxecc1=0.02 eccmean1=1e-4 eccscale1=eccmean1/3 #sigmask1=0.0003 #sigmask1=0.001

diskvisc1=1/sigmask1 beta1=0 ## ca -2.2 for outer planets ## if a is normal distribution meana1=1.0*math.sqrt(mstar1) scalea1=0.25*math.sqrt(mstar1) masspower1=3

numsteps=int((timemax-timebegin)/timestep)

snowline1=3*math.sqrt(mstar1)

frockdust=0.25

cdustrock=1/frockdust cgas=1/f

## switced off #taugas1=taugas1*1e6 #taudep1=taudep1*1e6 #tausolids1=taugas1

num1=int((aouter1-ainner1)/adelta1) disknum1=int((aouter1-ainner1)/diskdelta1) meanmass1=meanmassk1/num1 ## minimal planetesimal mass ! #aas1=np.linspace(ainner1, aouter1, num1) ## linear aas1=np.random.normal(loc=meana1, scale=scalea1, size=num1) ## normal distribution #aas1=np.random.random(num1)*(aouter1-ainner1)+ainner1 ## random #aas1=ainner1+np.power(np.linspace(0,(num1/5),num1 ), 1.2) #line1=np.linspace(0,num1,num1 ) #aas1=ainner1+0.4*np.power(1.15, line1 ) aas1=np.where(aas1<ainner1,ainner1*2, aas1) aas1=np.where(aas1>aouter1,aouter1/2, aas1)

len1=len(aas1)

snowyplanets1=np.copy(aas1) snowyplanets1=np.where(aas1>snowline1,cdustrock,1 )

adisk1=np.linspace(ainner1,aouter1, disknum1)

#dense1=np.copy(snowy1) #dense1=np.where(snowy1<1.1,1, 0.4)

#ms1=aas1*0.0+meanmass1 #ms1=np.random.rayleigh(meanmass1, len1) #vals1=np.random.rand(len1)*0+meanmassk1 ms1=np.random.normal(loc=meana1, scale=scalea1, size=num1) #ms1=np.power(vals1, 0)*meanmassk1 #ms1=np.power(vals1, masspower1)*meanmass1 ## ca 2-3 ? #ms1=0.0102*np.power(ms1, 2.7) ## < 1 % #ms1=np.exp(-ms1*500) #3 < 1 % #ms1=np.random.rayleigh(0.0033, num1)*meanmassk1 ## ca. 1 % bigger than 99% #ms1=vals1 ms1=np.where(ms1>maxmassk1,maxmassk1, ms1) ms1=np.where(ms1<minmassk1,minmassk1, ms1) ms1=ms1*snowyplanets1

mrocks1=np.copy(ms1) mices1=np.copy(ms1)*0 mgases1=np.copy(ms1)*0

#eccs1=1*abs(np.random.normal(loc=eccmean1, scale=eccscale1, size=len1))/ms1 eccs1=1*abs(np.random.normal(loc=eccmean1, scale=eccscale1, size=len1)) #omegas1=np.power(adisk1, -3/2) ## what is effect of mstar ?? #omegaplanetesimals1=np.power(aas1, -3/2) omegaplanetesimals1=np.sqrt(mstar1*np.power(aas1, -3)) vplanetesimals1=np.power(aas1, 2) ## nok denses1=np.copy(aas1) denses1=np.where(denses1<snowline1,1, 0.4) radiuses1=np.power((ms1/denses1), 0.27) vescapes1=np.sqrt(ms1/radiuses1) ## arbitrary unit! vrels1=vplanetesimals1*eccs1*eccs1 #sigmas1=np.power(aas1, beta1)*sigmask1

#quit(-1)

#collisionbetween=600 collisionbetween=10 #print("Initial values: ") #print (aas1) #print (ms1) #exit(-1)

time1=0 for ti in range(0, numsteps): len1=len(aas1) #if(len1<5): break omegas1=np.sqrt(mstar1*np.power(adisk1, -3)) ## sqrt(G*m/pow(a, 3)) snowy1=np.copy(adisk1) snowy1=np.where(adisk1>snowline1,cdustrock,1 ) beta1=0 #sigmagas1=np.power(adisk1, beta1)*sigmask1*cgas sigmagas1=np.random.normal(loc=meana1, scale=scalea1, size=num1)*sigmask1*cgas sigmagas1=sigmagas1*math.exp(-time1/taudep1) #sigmagas1=norm.pdf(adisk1, loc=0.9, scale =0.25)*sigmask1*cgas ## gas ring sigmasolidsmax1=sigmagas1*f*metallicity1 ## teor solids sigmadust1=sigmagas1*f*snowy1*metallicity1 ## fine grain anf pebble size dust or solids sigmarocks1=sigmagas1*f*frockdust*metallicity1 ## rocks sigmaices1=sigmagas1*f*snowy1*metallicity1-sigmarocks1 ## ices tdisk1=160*np.power(adisk1, -1/2) ## classic -6/7 vdisk1=np.sqrt(mstar1/adisk1) ## sqrt(G*m/r) vsoundh2=157.93*np.sqrt(tdisk1/2) ## molecular hydrogen vsoundhe2=157.93*np.sqrt(tdisk1/4)

if(ti==0): mplanetesimalsum1=np.sum(ms1) mdustsum1=np.sum(sigmadust1) #print("mpsum  : " , mplanetesimalsum1) #print("dustsum  : ",mdustsum1) #print("coalelescent fraction: ",(mplanetesimalsum1/(mplanetesimalsum1+mdustsum1)))

#quit(-1) #aas1, eccs1, ms1, sigmadust1, sigmagas1, sigmaices1, sigmarocks1=increase_mass_from_dustgas_1(mstar1, time1, timestep, taugas, tausolids, aas1, eccs1, ms1, adisk1, sigmadust1, sigmagas1, sigmaices1, sigmarocks1) aas1, eccs1, ms1, mrocks1, mices1, mgases1,sigmadust1, sigmagas1, sigmaices1, sigmarocks1=increase_mass_from_dustgas_1(mstar1,snowline1,incdustcoeff1, incgascoeff1, bhill1, time1, timestep, taugas1, tausolids1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, adisk1, sigmadust1, sigmagas1, sigmaices1, sigmarocks1)

if(np.any(ms1)>6000): print("Mass over limit") break

#print(time1, ms1) len1=len(aas1) #if(len1<5): break

if(len1>maxplanets): if((ti%collisionbetween)==0): idx1=np.random.randint(len1) # #print("*") #aas1, eccs1, ms1, mrocks1, mices1, mgases1, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals(mstar1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, idx1, bee1, adisk1, sigmadust1, sigmaices1, sigmarocks1) aas1, eccs1, ms1, mrocks1, mices1, mgases1, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals_from_ecc_and_hill_radius(1, mstar1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, idx1, bee1, adisk1, sigmadust1, sigmaices1, sigmarocks1)

eccs1=eccs1*1.5 if(np.any(eccs1>0.1)): eccs1=eccs1/2 collisionbetween=int(collisionbetween*1.01) # #print(aas1) # #print(ms1)

aas1=aas1*migk1

#eccs1=eccs1*0+

time1=time1+timestep

len1=len(aas1) if(len1<5): mummu=1

# print("---") # print(time1) # print(aas1) # print(eccs1) # print(ms1) len1=len(aas1) #if(len1<5): break

aas3=aas1 #print(".")

eccs3=eccs1 ms3=ms1 mrocks3=mrocks1 mices3=mices1 mgases3=mgases1

#print("Coalesced0: ") #print (aas3) #print(eccs3) #print(ms3) #plt.scatter(aas3, aas3*0+0.008, s=ms3*2000, color="g", facecolor="none", lw=5, alpha=0.5)

#aas3, eccs3, ms3, sigmadust1, sigmaices1, sigmarocks1=combine_small_planetesimals_to_bigger_objects(aas3, eccs3, ms3, mlimit1, adisk1, sigmadust1, sigmaices1, sigmarocks1) aas3, eccs3, ms3, mrocks3, mices3, mgases3,sigmadust1, sigmaices1, sigmarocks1=combine_small_planetesimals_to_bigger_objects(aas3, eccs3, ms3,mrocks3, mices3, mgases3, mlimit1, adisk1, sigmadust1, sigmaices1, sigmarocks1)

#print("Coalesced1: ") #print (aas3) #print(eccs3) #print(ms3)

#print("Coalesced total1 ", np.sum(ms3), " ME to ",len(ms3) ," planets" )

#exit(-1)

relativedistance1=0.45 ## ## nok !!! aas3, eccs3, ms3, mrocks3, mices3, mgases3=combine_by_relative_distance_1(aas3, eccs3, ms3,mrocks3, mices3, mgases3,relativedistance1) #plt.scatter(aas3, aas3*0+0.004, s=ms3*1000, color="k", facecolor="r", alpha=0.5)

print("Coalesced3: ") print (aas3) #print(eccs3) print(ms3) print("Coalesced total3 ", np.sum(ms3), " ME to ",len(ms3) ," planets" )

#plt.show()

#plt.plot(aas2, ms2) #plt.plot(aas1, sigmas1) return(aas3,eccs3, ms3, mrocks3, mices3, mgases3)

def generate_outer_planets_1(mstar1, surf1, metallicity1, alpha1, taugas1, taudep1, lifetime1, migratek1): ## gas giant maxplanets=8

#sigmask1=0.1*surf1 bhill1=2.4 ## 3.5 by default add to more planets bee1=3.5 f=0.013*metallicity1 frockdust=0.25

snowline1=3*math.sqrt(mstar1) sigmask1=0.001*surf1 ## dust and gas add coefficiens incdustcoeff1=1.063 ## 0.22*4 incgascoeff1=900

ainner1=5*math.sqrt(mstar1) aouter1=50*math.sqrt(mstar1)

#alpha1=0.001

migk1=migratek1 ## migration 1 no migration, 0 falls to star!

taugas1=taugas1*1e6 taudep1=taudep1*1e6

tausolids1=taugas1 ## not used

timebegin=0 timemax=int(lifetime1*1e6) timestep=1000

#sigmask1=0.02*12 #sigmask1=0.02

mlimit1=0.02 ## min planet mass #bee1=3.5

adelta1=0.3*math.sqrt(mstar1) diskdelta1=0.3*math.sqrt(mstar1) #meanmassk1=2e-7*2 meanmassk1=1e-8*4 ## stonyicy mean mass, stony planetesimals, icy has ca 4x maxmassk1=meanmassk1*1000 minmassk1=meanmassk1*1 maxecc1=0.03 eccmean1=1e-4 eccscale1=eccmean1/3 #sigmask1=0.0003 #sigmask1=0.001

diskvisc1=1/sigmask1 beta1=0 ## ca -2.2 for outer planets ## if a is normal distribution meana1=1*math.sqrt(mstar1) scalea1=0.25*math.sqrt(mstar1) masspower1=3

numsteps=int((timemax-timebegin)/timestep)

cdustrock=1/frockdust cgas=1/f

num1=int((aouter1-ainner1)/adelta1) disknum1=int((aouter1-ainner1)/diskdelta1) meanmass1=meanmassk1/num1 ## minimal planetesimal mass ! aas1=np.linspace(ainner1, aouter1, num1) ## linear #aas1=np.random.normal(loc=meana1, scale=scalea1, size=num1) ## normal distribution #aas1=np.random.random(num1)*(aouter1-ainner1)+ainner1 ## random #aas1=ainner1+np.power(np.linspace(0,(num1/5),num1 ), 1.2) #line1=np.linspace(0,num1,num1 ) #aas1=ainner1+0.4*np.power(1.15, line1 ) aas1=np.where(aas1<ainner1,ainner1*2, aas1) aas1=np.where(aas1>aouter1,aouter1/2, aas1)

len1=len(aas1)

snowyplanets1=np.copy(aas1) snowyplanets1=np.where(aas1>snowline1,cdustrock,1 )

adisk1=np.linspace(ainner1,aouter1, disknum1)

#dense1=np.copy(snowy1) #dense1=np.where(snowy1<1.1,1, 0.4)

#ms1=aas1*0.0+meanmass1 #ms1=np.random.rayleigh(meanmass1, len1) vals1=np.random.rand(len1)*0+meanmassk1 #ms1=np.random.normal(loc=meana1, scale=scalea1, size=num1) #ms1=np.power(vals1, 0)*meanmassk1 #ms1=np.power(vals1, masspower1)*meanmass1 ## ca 2-3 ? #ms1=0.0102*np.power(ms1, 2.7) ## < 1 % #ms1=np.exp(-ms1*500) #3 < 1 % #ms1=np.random.rayleigh(0.0033, num1)*meanmassk1 ## ca. 1 % bigger than 99% ms1=np.power(vals1, 0)*meanmassk1 ms1=np.where(ms1>maxmassk1,maxmassk1, ms1) ms1=np.where(ms1<minmassk1,minmassk1, ms1) ms1=ms1*snowyplanets1

mrocks1=np.copy(ms1) mices1=np.copy(ms1)*0 mgases1=np.copy(ms1)*0 #eccs1=1*abs(np.random.normal(loc=eccmean1, scale=eccscale1, size=len1))/ms1 eccs1=1*abs(np.random.normal(loc=eccmean1, scale=eccscale1, size=len1)) #omegas1=np.power(adisk1, -3/2) ## what is effect of mstar ?? #omegaplanetesimals1=np.power(aas1, -3/2) omegaplanetesimals1=np.sqrt(mstar1*np.power(aas1, -3)) vplanetesimals1=np.power(aas1, 2) ## nok denses1=np.copy(aas1) denses1=np.where(denses1<snowline1,1, 0.4) radiuses1=np.power((ms1/denses1), 0.27) vescapes1=np.sqrt(ms1/radiuses1) ## arbitrary unit! vrels1=vplanetesimals1*eccs1*eccs1 #sigmas1=np.power(aas1, beta1)*sigmask1

#quit(-1)

#collisionbetween=100 # 600 collisionbetween=8 #print("Initial values: ") #print (aas1) #print (ms1) #exit(-1)

time1=0 for ti in range(0, numsteps): len1=len(aas1) #if(len1<5): break omegas1=np.sqrt(mstar1*np.power(adisk1, -3)) ## sqrt(G*m/pow(a, 3)) snowy1=np.copy(adisk1) snowy1=np.where(adisk1>snowline1,cdustrock,1 ) #beta1=-2.2 #beta1=-2.5 beta1=-1 #sigmagas1=np.power(adisk1, beta1)*sigmask1*cgas #sigmagas1=np.copy(adisk1)*0+1*sigmask1*cgas #sigmagas1=np.random.normal(loc=meana1, scale=scalea1, size=num1)*sigmask1*cgas #sigmagas1=sigmagas1*math.exp(-time1/taugas) #sigmagas1=norm.pdf(adisk1, loc=0.9, scale =0.25)*sigmask1*cgas ## gas ring macc1=calc_macc_1(mstar1, taugas1, time1) sigmagas1=calculate_sigmas_from_macc_0(mstar1, macc1, adisk1, alpha1) sigmagas1=sigmagas1*math.exp(-time1/taudep1)

sigmasolidsmax1=sigmagas1*f*metallicity1 ## teor solids sigmadust1=sigmagas1*f*snowy1*metallicity1 ## fine grain anf pebble size dust or solids sigmarocks1=sigmagas1*f*frockdust*metallicity1 ## rocks sigmaices1=sigmagas1*f*snowy1*metallicity1-sigmarocks1 ## ices tdisk1=160*np.power(adisk1, -1/2) ## classic -6/7 vdisk1=np.sqrt(mstar1/adisk1) ## sqrt(G*m/r) vsoundh2=157.93*np.sqrt(tdisk1/2) ## molecular hydrogen vsoundhe2=157.93*np.sqrt(tdisk1/4)

if(ti==0): mplanetesimalsum1=np.sum(ms1) mdustsum1=np.sum(sigmadust1) #print("mpsum  : " , mplanetesimalsum1) #print("dustsum  : ",mdustsum1) #print("coalelescent fraction: ",(mplanetesimalsum1/(mplanetesimalsum1+mdustsum1)))

#quit(-1) aas1, eccs1, ms1, mrocks1, mices1, mgases1,sigmadust1, sigmagas1, sigmaices1, sigmarocks1=increase_mass_from_dustgas_1(mstar1, snowline1, incdustcoeff1, incgascoeff1, bhill1, time1, timestep, taugas1, tausolids1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, adisk1, sigmadust1, sigmagas1, sigmaices1, sigmarocks1)

if(np.any(ms1)>6000): print("Mass over limit") break

#print(time1, ms1) len1=len(aas1) #if(len1<5): break

if(len1>maxplanets): if((ti%collisionbetween)==0): idx1=np.random.randint(len1) #print("*") #aas1, eccs1, ms1, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals(mstar1, aas1, eccs1, ms1, idx1, bee1, adisk1, sigmadust1, sigmaices1, sigmarocks1) #aas1, eccs1, ms1, mrocks1, mices1, mgases1, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals(mstar1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, idx1, bee1, adisk1, sigmadust1, sigmaices1, sigmarocks1) aas1, eccs1, ms1, mrocks1, mices1, mgases1, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals_from_ecc_and_hill_radius(1,mstar1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, idx1, bee1, adisk1, sigmadust1, sigmaices1, sigmarocks1)

eccs1=eccs1*1.000 if(np.any(eccs1>0.05)): eccs1=eccs1/2 collisionbetween=int(collisionbetween*1.01) # #print(aas1) # #print(ms1)

aas1=aas1*migk1

#eccs1=eccs1*0+

time1=time1+timestep

len1=len(aas1) if(len1<10): mummu=1 #print("---") #print(time1) #print(aas1) #print(eccs1) #print(ms1) len1=len(aas1) #if(len1<5): break

aas3=aas1 eccs3=eccs1 ms3=ms1 mrocks3=mrocks1 mices3=mices1 mgases3=mgases1

#print("Coalesced0: ") #print (aas3) #print(eccs3) #print(ms3) #print(mrocks3) #print(mices3) #print(mgases3)

#quit(-1)

#plt.scatter(aas3, aas3*0+0.008, s=ms3*2000, color="g", facecolor="none", lw=5, alpha=0.5)

aas3, eccs3, ms3,mrocks3, mices3, mgases3,sigmadust1, sigmaices1, sigmarocks1=combine_small_planetesimals_to_bigger_objects(aas3, eccs3, ms3,mrocks3, mices3, mgases3, mlimit1, adisk1, sigmadust1, sigmaices1, sigmarocks1) #print("Coalesced1: ") #print (aas3) #print(eccs3) #print(ms3) #print("rocks, ices, gases: ") #print(mrocks3) #print(mices3) #print(mgases3)

#print (len(aas3), len(eccs3), len(ms3), len(mrocks3), len(mices3), len(mgases3))

#quit(-1)

#print("Coalesced total1 ", np.sum(ms3), " ME to ",len(ms3) ," planets" )

#exit(-1)

#print("Reldist ...") relativedistance1=0.5 ## ## nok !!! aas3, eccs3, ms3, mrocks3, mices3, mgases3 =combine_by_relative_distance_1(aas3, eccs3, ms3,mrocks3, mices3, mgases3, relativedistance1)

#print(" aas ms") #print(aas3) #print(ms3) #print("rocks, ices, gases: ") #print(mrocks3) #print(mices3) #print(mgases3)

#print (len(aas3), len(eccs3), len(ms3), len(mrocks3), len(mices3), len(mgases3)) #plt.scatter(aas3, aas3*0+0.004, s=ms3*1000, color="k", facecolor="r", alpha=0.5) #plt.show() #quit(-1)

#plt.scatter(aas3, aas3*0+0.004, s=ms3*1000, color="k", facecolor="r", alpha=0.5) #aext1=5.2 ## NOK #aas3=find_nearest_resonances_by_external_1(aas3, ms3, aext1)

#print("Coalesced3: ") #print (aas3) #print(eccs3) #print(ms3) #print("Coalesced total3 ", np.sum(ms3), " ME to ",len(ms3) ," planets" )

#plt.show()

#plt.plot(aas2, ms2) #plt.plot(aas1, sigmas1) return(aas3,eccs3, ms3, mrocks3, mices3, mgases3)

def generate_planet_system_type_1(seed1, mstar1, surf1, metallicity1, alpha1, taugas1, taudep1, lifetime1, migratek1): np.random.seed(seed1)

## generate outer planets aas2, eccs2, ms2, mrocks2, mices2, mgases2=generate_outer_planets_1(mstar1, surf1, metallicity1, alpha1, taugas1, taudep1, lifetime1, migratek1) #print (len(aas2), len(eccs2), len(ms2), len(mrocks2), len(mices2), len(mgases2))

print(aas2) print(ms2) #print(mrocks2) #print(mices2) #print(mgases2) #quit(-1)

## generate inner planets

lifetime2=50 aas1, eccs1, ms1, mrocks1, mices1, mgases1=generate_inner_planets_1(mstar1, surf1, metallicity1, taugas1, taudep1, lifetime2, migratek1)

## attempt to adjust resonances to outer planets ...

jdex2=np.where(ms2>10)[0][0]

aext2=aas2[jdex2] #aext2=5.2

## NOK aas1=find_nearest_resonances_by_external_1(aas1, ms1, aext2) print("Inner system 0: ") print(aas1) print(ms1) #print(mrocks1) #print(mices1) #print(mgases1) #quit(-1)

#print (len(aas1), len(eccs1), len(ms1), len(mrocks1), len(mices1), len(mgases1)) #quit(-1) aas3=np.append(aas1, aas2) eccs3=np.append(eccs1, eccs2) ms3=np.append(ms1, ms2) mrocks3=np.append(mrocks1, mrocks2) mices3=np.append(mices1, mices2) mgases3=np.append(mgases1, mgases2)

## combine giant planets that are massive, 3.5 hill radius

len1=len(aas3)

for n in range(0,100): idx1=np.random.randint(len1) beibe1=3.5 aas3, eccs3, ms3, mrocks3, mices3, mgases3, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals_from_ecc_and_hill_radius(1,mstar1, aas3, eccs3, ms3, mrocks3, mices3, mgases3, idx1, beibe1, None, None, None, None)

aas3, eccs3, ms3, mrocks3, mices3, mgases3=process_inner_planet_system_by_stability(mstar1, aas3, eccs3, ms3, mrocks3, mices3, mgases3)

#plotting3(aas1, eccs1,ms1, 0, 3) #plotting3(aas2, eccs2,ms2, 5, 40)

#plt.show()

return(aas3, eccs3, ms3, mrocks3, mices3, mgases3)

def generate_super_planets_1(mstar1, surf1, metallicity1, alpha1, taugas1, taudep1, lifetime1, migratek1):

print ("super planets ...")

bhill1=20 bee1=20

## gas giant maxplanets=8

sigmask1=0.0001*surf1 ## dust and gas add coefficiens incdustcoeff1=0.22 incgascoeff1=900

f=0.013*metallicity1 frockdust=0.25

#alpha1=0.001

migk1=migratek1 ## migration 1 no migration, 0 falls to star! taugas1=taugas1*1e6 taudep1=taudep1*1e6

tausolids1=taugas1 ## not used timebegin=0 timemax=int(lifetime1*1e6) timestep=1000

snowline1=3*math.sqrt(mstar1)

#sigmask1=0.02*12 #sigmask1=0.02

mlimit1=0.02 ## min planet mass maxecc1=0.001

ainner1=5*math.sqrt(mstar1) aouter1=200*math.sqrt(mstar1) adelta1=0.1*math.sqrt(mstar1) diskdelta1=0.1*math.sqrt(mstar1) #meanmassk1=2e-7*2 meanmassk1=1e-7*4 ## stonyicy mean mass, stony planetesimals, icy has ca 4x maxmassk1=meanmassk1*1000 minmassk1=meanmassk1*1

eccmean1=1e-4 eccscale1=eccmean1/3 #sigmask1=0.0003 #sigmask1=0.001

diskvisc1=1/sigmask1 beta1=0 ## ca -2.2 for outer planets ## if a is normal distribution meana1=1*math.sqrt(mstar1) scalea1=0.25*math.sqrt(mstar1) masspower1=3

numsteps=int((timemax-timebegin)/timestep)

cdustrock=1/frockdust cgas=1/f

num1=int((aouter1-ainner1)/adelta1) disknum1=int((aouter1-ainner1)/diskdelta1) meanmass1=meanmassk1/num1 ## minimal planetesimal mass ! aas1=np.linspace(ainner1, aouter1, num1) ## linear #aas1=np.random.normal(loc=meana1, scale=scalea1, size=num1) ## normal distribution #aas1=np.random.random(num1)*(aouter1-ainner1)+ainner1 ## random #aas1=ainner1+np.power(np.linspace(0,(num1/5),num1 ), 1.2) #line1=np.linspace(0,num1,num1 ) #aas1=ainner1+0.4*np.power(1.15, line1 ) aas1=np.where(aas1<ainner1,ainner1*2, aas1) aas1=np.where(aas1>aouter1,aouter1/2, aas1)

len1=len(aas1)

snowyplanets1=np.copy(aas1) snowyplanets1=np.where(aas1>snowline1,cdustrock,1 )

adisk1=np.linspace(ainner1,aouter1, disknum1)

#dense1=np.copy(snowy1) #dense1=np.where(snowy1<1.1,1, 0.4)

#ms1=aas1*0.0+meanmass1 #ms1=np.random.rayleigh(meanmass1, len1) vals1=np.random.rand(len1)*0+meanmassk1 #ms1=np.random.normal(loc=meana1, scale=scalea1, size=num1) #ms1=np.power(vals1, 0)*meanmassk1 #ms1=np.power(vals1, masspower1)*meanmass1 ## ca 2-3 ? #ms1=0.0102*np.power(ms1, 2.7) ## < 1 % #ms1=np.exp(-ms1*500) #3 < 1 % #ms1=np.random.rayleigh(0.0033, num1)*meanmassk1 ## ca. 1 % bigger than 99% ms1=np.power(vals1, 0)*meanmassk1 ms1=np.where(ms1>maxmassk1,maxmassk1, ms1) ms1=np.where(ms1<minmassk1,minmassk1, ms1) ms1=ms1*snowyplanets1

mrocks1=np.copy(ms1) mices1=np.copy(ms1)*0 mgases1=np.copy(ms1)*0 #eccs1=1*abs(np.random.normal(loc=eccmean1, scale=eccscale1, size=len1))/ms1 eccs1=1*abs(np.random.normal(loc=eccmean1, scale=eccscale1, size=len1)) #omegas1=np.power(adisk1, -3/2) ## what is effect of mstar ?? #omegaplanetesimals1=np.power(aas1, -3/2) omegaplanetesimals1=np.sqrt(mstar1*np.power(aas1, -3)) vplanetesimals1=np.power(aas1, 2) ## nok denses1=np.copy(aas1) denses1=np.where(denses1<snowline1,1, 0.4) radiuses1=np.power((ms1/denses1), 0.27) vescapes1=np.sqrt(ms1/radiuses1) ## arbitrary unit! vrels1=vplanetesimals1*eccs1*eccs1 #sigmas1=np.power(aas1, beta1)*sigmask1

#quit(-1)

#collisionbetween=100 # 600 collisionbetween=1 ## 10 #print("Initial values: ") #print (aas1) #print (ms1) #exit(-1)

time1=0 for ti in range(0, numsteps): len1=len(aas1) #if(len1<5): break omegas1=np.sqrt(mstar1*np.power(adisk1, -3)) ## sqrt(G*m/pow(a, 3)) snowy1=np.copy(adisk1) snowy1=np.where(adisk1>snowline1,cdustrock,1 ) #beta1=-2.2 #beta1=-2.5 beta1=0 #sigmagas1=np.power(adisk1, beta1)*sigmask1*cgas #sigmagas1=np.copy(adisk1)*0+1*sigmask1*cgas #sigmagas1=np.random.normal(loc=meana1, scale=scalea1, size=num1)*sigmask1*cgas

#sigmagas1=norm.pdf(adisk1, loc=0.9, scale =0.25)*sigmask1*cgas ## gas ring macc1=calc_macc_1(mstar1, taugas1, time1) sigmagas1=calculate_sigmas_from_macc_0(mstar1, macc1, adisk1, alpha1) sigmagas1=sigmagas1*math.exp(-time1/taudep1)

sigmasolidsmax1=sigmagas1*f*metallicity1 ## teor solids sigmadust1=sigmagas1*f*snowy1*metallicity1 ## fine grain anf pebble size dust or solids sigmarocks1=sigmagas1*f*frockdust*metallicity1 ## rocks sigmaices1=sigmagas1*f*snowy1*metallicity1-sigmarocks1 ## ices tdisk1=160*np.power(adisk1, -1/2) ## classic -6/7 vdisk1=np.sqrt(mstar1/adisk1) ## sqrt(G*m/r) vsoundh2=157.93*np.sqrt(tdisk1/2) ## molecular hydrogen vsoundhe2=157.93*np.sqrt(tdisk1/4)

if(ti==0): mplanetesimalsum1=np.sum(ms1) mdustsum1=np.sum(sigmadust1) print("mpsum  : " , mplanetesimalsum1) print("dustsum  : ",mdustsum1) print("coalescent fraction: ",(mplanetesimalsum1/(mplanetesimalsum1+mdustsum1)))

#quit(-1) #aas1, eccs1, ms1, mrocks1, mices1, mgases1,sigmadust1, sigmagas1, sigmaices1, sigmarocks1=increase_mass_from_dustgas_1(1,mstar1, snowline1, incdustcoeff1, incgascoeff1, bhill1, time1, timestep, taugas1, tausolids1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, adisk1, sigmadust1, sigmagas1, sigmaices1, sigmarocks1) aas1, eccs1, ms1, mrocks1, mices1, mgases1,sigmadust1, sigmagas1, sigmaices1, sigmarocks1=increase_mass_from_dustgas_1(mstar1,snowline1, incdustcoeff1, incgascoeff1, bhill1, time1, timestep, taugas1, tausolids1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, adisk1, sigmadust1, sigmagas1, sigmaices1, sigmarocks1)

if(np.any(ms1)>6000): print("Mass over limit") break

#print(time1, ms1) len1=len(aas1) #if(len1<5): break

if(len1>maxplanets): if((ti%collisionbetween)==0): idx1=np.random.randint(len1) #print("*") #aas1, eccs1, ms1, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals(mstar1, aas1, eccs1, ms1, idx1, bee1, adisk1, sigmadust1, sigmaices1, sigmarocks1) #aas1, eccs1, ms1, mrocks1, mices1, mgases1, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals(mstar1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, idx1, bee1, adisk1, sigmadust1, sigmaices1, sigmarocks1) aas1, eccs1, ms1, mrocks1, mices1, mgases1, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals_from_ecc_and_hill_radius(1,mstar1, aas1, eccs1, ms1, mrocks1, mices1, mgases1, idx1, bee1, adisk1, sigmadust1, sigmaices1, sigmarocks1)

eccs1=eccs1*1.000 if(np.any(eccs1>0.05)): eccs1=eccs1/2 collisionbetween=int(collisionbetween*1.01) # #print(aas1) # #print(ms1)

aas1=aas1*migk1

#eccs1=eccs1*0+

time1=time1+timestep

len1=len(aas1) if(len1<5): mummu=1 #print("---") #print(time1) #print(aas1) #print(eccs1) #print(ms1) len1=len(aas1) #if(len1<5): break

aas3=aas1 eccs3=eccs1 ms3=ms1 mrocks3=mrocks1 mices3=mices1 mgases3=mgases1

print("Coalesced0: ") print (aas3) print(eccs3) print(ms3) print(mrocks3) print(mices3) print(mgases3)

#quit(-1)

#plt.scatter(aas3, aas3*0+0.008, s=ms3*2000, color="g", facecolor="none", lw=5, alpha=0.5)

aas3, eccs3, ms3,mrocks3, mices3, mgases3,sigmadust1, sigmaices1, sigmarocks1=combine_small_planetesimals_to_bigger_objects(aas3, eccs3, ms3,mrocks3, mices3, mgases3, mlimit1, adisk1, sigmadust1, sigmaices1, sigmarocks1) print("Coalesced1: ") print (aas3) print(eccs3) print(ms3) print("rocks, ices, gases: ") print(mrocks3) print(mices3) print(mgases3)

print (len(aas3), len(eccs3), len(ms3), len(mrocks3), len(mices3), len(mgases3))

#quit(-1)

#print("Coalesced total1 ", np.sum(ms3), " ME to ",len(ms3) ," planets" )

#exit(-1)

print("Reldist ...") relativedistance1=0.5 ## ## nok !!! #"" sloow aas3, eccs3, ms3, mrocks3, mices3, mgases3 =combine_by_relative_distance_1(aas3, eccs3, ms3,mrocks3, mices3, mgases3, relativedistance1)

print(" aas ms") print(aas3) print(ms3) print("rocks, ices, gases: ") print(mrocks3) print(mices3) print(mgases3)

len1=len(aas3) ## combine giats etc.... for n in range(0,100): idx1=np.random.randint(len1) beibe1=3.5 aas3, eccs3, ms3, mrocks3, mices3, mgases3, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals_from_ecc_and_hill_radius(1,mstar1, aas3, eccs3, ms3, mrocks3, mices3, mgases3, idx1, beibe1, None, None, None, None)

print (len(aas3), len(eccs3), len(ms3), len(mrocks3), len(mices3), len(mgases3)) #plt.scatter(aas3, aas3*0+0.004, s=ms3*1000, color="k", facecolor="r", alpha=0.5) #plt.show() #quit(-1)

#plt.scatter(aas3, aas3*0+0.004, s=ms3*1000, color="k", facecolor="r", alpha=0.5) #aext1=5.2 ## NOK #aas3=find_nearest_resonances_by_external_1(aas3, ms3, aext1)

print("Coalesced3: ") print (aas3) print(eccs3) print(ms3) print("Coalesced total3 ", np.sum(ms3), " ME to ",len(ms3) ," planets" )

#plt.show()

#plt.plot(aas2, ms2) #plt.plot(aas1, sigmas1) return(aas3,eccs3, ms3, mrocks3, mices3, mgases3)

def generate_planet_system_type_2(seed1, mstar1, surf1, metallicity1, alpha1, taugas1,taudep1, lifetime1, migratek1):

np.random.seed(seed1)

print ("Generating planet system type 2 ...") ## generate outer planets

aas3, eccs3, ms3, mrocks3, mices3, mgases3=generate_super_planets_1(mstar1, surf1, metallicity1, alpha1, taugas1, taudep1, lifetime1, migratek1)

sigmadust1=None sigmaices1=None sigmarocks1=None

len1=len(aas3) ## combine giats etc.... for n in range(0,100): idx1=np.random.randint(len1) beibe1=3.5 aas3, eccs3, ms3, mrocks3, mices3, mgases3, sigmadust1, sigmaices1, sigmarocks1=collect_nearby_planetesimals_from_ecc_and_hill_radius(1,mstar1, aas3, eccs3, ms3, mrocks3, mices3, mgases3, idx1, beibe1, None, None, None, None)

return(aas3, eccs3, ms3, mrocks3, mices3, mgases3)

1. 2. main program ...

1. 1. TODO

1. 1. eccentricity pump-up for small and very big planets
   2. then combine if necessary

1. 1. more accurate temperatures of planets

1. 1. diskmass?
   2. disksize?
   3. taudust1??
   4. disk h/r
   5. gasdensity, giant planet blocks

1. 1. omega effect of star mss
   2. sound speed
   3. pebble accretion ?
   4. resonances

1. seed1, mstar1=read_command_line()

1. 1. basic control parameters of planet system generator test

1. seed1=None ## random system

seed1=88 #3 three planets

1. seed1=777

1. seed1=999999

1. seed1=123123

1. seed1=54321 ## ok

1. seed1=4444 ## many planets

1. seed1=3

1. seed1=None

1. seed1=9999 ##

1. 1. par1=sys.argv[1]

1. random_status0= np.random.get_state()

1. gotseed1=np.random.get_state()[1][1]

1. print(seed1, gotseed1)

1. print(random_status0)

1. quit(-1)

1. print(tau1, teff1, tsurf1, tsurf1-273.15)

1. quit(-1)

name1="Planets of Wasa" outfilename1='planets_0000.csv'

planet_system_type1=1 ## 1 two sets like solar system rocky and icy/gas 2: some wet/icy/gassy, like Trappist or so on

drawskala1=1 ## 1 normal drawing scala of planets near their star.

mstar1=1 ## 1 sun surf1=1.0 ## 1 normal metallicity1=1 ## 1 sun alpha1=1.0e-3 ## tex 1e-2, 1e-3, 1e-4 viscosity alpha for giant planets only! taugas1=6 ## #3 viscosity tau of disk ca 15 normal 1 small planets taudep1=3 ## 15 depletion tau of disk lifetime1=15 ## ca 15 normal 1 small planets

migratek1=0.0 ## 0 default migration: 0 no migration

1. 1. gaseous and rocky, inner and outer planets

1. 1. read command line ...
   2. tex python planetsysgen1 -s 12345
   3. uncomment line below in you want command line arguments

1. name1, outfilename1, type1, seed1, mstar1, metallicity1, surfacedensity1,lifetime1, taugas1, taudep1, alpha1, migratek1=read_command_line()

1. 1. to a coefficient

migratek2=1-migratek1

if (planet_system_type1==0): planet_system_type1=1

if (planet_system_type1==1): print("Attempting to make two sets of planets: inner rocky and outer icy/gaseous planets.") aas1, eccs1, ms1, mrocks1, mices1, mgases1=generate_planet_system_type_1(seed1, mstar1, surf1, metallicity1,alpha1, taugas1, taudep1, lifetime1, migratek2)

if(planet_system_type1==2): print("Making one set icy/wet/gassy planets.") ## "similar" planets, one far formed set only ## so icy or wet super-earths , neptunes. gaseous ... aas1, eccs1, ms1, mrocks1, mices1, mgases1=generate_planet_system_type_2(seed1, mstar1, surf1, metallicity1,alpha1, taugas1, taudep1, lifetime1, migratek2)

1. lstar1=math.pow(mstar1, 3.5) ## very coarse est

lstar1=calc_lstar_from_mass(mstar1)

teffstar1=5778*math.pow(mstar1,0.54)

1. time_main_sequence1_gyr1=10.9*mstar1/lstar1
2. time_main_sequence1=math.pow()
   1. earsth correspondingh habitable zone

ahz1=math.sqrt(lstar1)*math.pow((teffstar1/5780), -0.799)

mall1=mrocks1+mices1+mgases1

rockp1=mrocks1/mall1 icep1=mices1/mall1 gasp1=mgases1/mall1

rs1=np.copy(aas1)*0

for n in range(0, len(ms1)): rs1[n]=radius_s1(ms1[n], rockp1[n], icep1[n], gasp1[n])

co2pressure=400e-6 h2opressure=0 solarconstant=2*1368 albedo=0.3

tau1=calc_greenhouse_tau_co2_h2o(co2pressure, h2opressure)

teff1=calcu_temperature(solarconstant, albedo,1) tsurf1=teff1*(1+3/4*tau1)

rotperiod=24 orbitperiods1=np.sqrt((aas1*aas1*aas1)/ms1 ) temps1=planet_temperatures_1(mstar1, lstar1, aas1) correspondic_solar_distances1=aas1/ahz1 sols1=lstar1/np.power(aas1,2)

minimol1=minimum_molecule_mass_to_hold_in_atmosphere(temps1, ms1, rs1) locked1=is_tidal_locked(5,mstar1,ms1, aas1, 12) patms1=atmosphere_pressure_estimation(ms1, rs1, temps1)

1. 1. attempt to calculate greenhouse effect

1. 1. from hat greenhouse effect estimation

temps2=temps1*1.126*np.log10(patms1) for n in range(0, len(temps1)): if(ms1[n]<20): if(temps2[n]>temps1[n]): temps1[n]=temps2[n]

1. 1. maybe better ....

co2pressures1=400e-6*patms1 h2opressures1=0*patms1 solarconstants1=sols1*1368 albedos1=0.00

tauplanets1=calc_greenhouse_tau_co2_h2o(co2pressures1, h2opressures1)

teffplanets1=calcu_temperature(solarconstants1, albedos1,1) tsurfplanets1=teffplanets1*(1+3/4*tauplanets1)

len1=len(tsurfplanets1)

1. 1. venusian like worlds ...

for n in range(0, len1): t1=tsurfplanets1[n] p1=patms1[n] m1=ms1[n] if(ms1[n]>0.6): if(t1>373): ## runaway greenhouse patms1[n]=patms1[n]*m1*80 co2pressures1=patms1[n]*0.96

tauplanets1=calc_greenhouse_tau_co2_h2o(co2pressures1, h2opressures1) tsurfplanets1=teffplanets1*(1+3/4*tauplanets1)

esis1=calc_esis(ms1, rs1, tsurfplanets1)

print(temps1)

1. print(temps2)

print(tsurfplanets1)

1. quit(-1)

gees1=gravities_earths(ms1, rs1) vescs1=vescs_earths(ms1, rs1)*11.86

rots1=rotation_periods_hours(ms1, rs1)

magnetic1=magnetic_potential_estimation(ms1, rs1, rots1)

planettypes1=np.copy(aas1)*0

len1=len(aas1)

for n in range(0, len(planettypes1)): #rs1[n]=radius_s1(ms1[n], rockp1[n], icep1[n], gasp1[n]) planettypes1[n]=10 ## rocky by default ## coarse estimation of stellar flux #sol1=1/math.pow(aas1[n],2)*math.sqrt(math.pow(mstar1,4)) ## estim check this ## is it Terra-like #if(sol1>0.85): # if(sol1<1.25): # if(ms1[n]>0.5): # if(ms1[n]<2): # planettypes1[n]=6 if(esis1[n]>0.97): planettypes1[n]=12

if(icep1[n]>0.02): planettypes1[n]=2 ## icy if(ms1[n]>12): planettypes1[n]=3 ## neptune if(ms1[n]>30): planettypes1[n]=4 ## saturn if(ms1[n]>150): planettypes1[n]=5 ## jupiter

if (rockp1[n]>99.9): if(tsurfplanets1[n]>450): if(ms1[n]>0.6): planettypes1[n]=11 ## venusian

planettypes1[n]=1 ## rocky desert by default if(ms1[n]<0.07): planettypes1[n]=10 ## stony or mercurial

if(tsurfplanets1[n]>340): ## warm terran if(ms1[n]>0.6): planettypes1[n]=14

if(sols1[n]<0.8): planettypes1[n]=13 ## martian

if(ms1[n]>0.06): ## terran size

if(sols1[n]>1.05): if(sols1[n]<1.36): planettypes1[n]=14## warm terran

if(esis1[n]>0.60): planettypes1[n]=15 ## near terran

if(esis1[n]>0.94): planettypes1[n]=12 ## terran

if(sols1[n]>1.3): planettypes1[n]=11 ## venusian

if (icep1[n]>0.02): if(gasp1[n]<0.01): planettypes1[n]=2 ## icy by default if(ms1[n]<0.07): planettypes1[n]=10 ## stony or mercurial if(sols1[n]>0.95): planettypes1[n]=21 ## waterworld if(sols1[n]>1.3): planettypes1[n]=21 ## venusian waterworld if(sols1[n]<0.8): planettypes1[n]=22 #icy

if(gasp1[n]>0.01): if(ms1[n]>0.1): planettypes1[n]=32 ## gasdwarf or mini-neptune

if(ms1[n]>13): planettypes1[n]=3 ## neptune if(ms1[n]>30): planettypes1[n]=4 ## saturn if(ms1[n]>150): planettypes1[n]=5 ## jupiter

len1=len(rots1)

for n in range(0, len1): if (locked1[n]==1): rots1[n]=orbitperiods1[n]*24*365

1. quit(-1)

1. rockp1=mrocks1/ms1
2. icep1=mices1/ms1
3. gasp1=mgases1/ms1

1. print(ms1)
2. print(mall1)
3. print(rockp1)
4. print(icep1)
5. print(gasp1)

print() print() print("Main input parameters") print("-----------------")

1. print("seed ",random_status0)

print("planet system name ", name1) print("planet system type ", planet_system_type1) print("mstar ", mstar1) print("surface density ",surf1) print("metallicity ",metallicity1) print("disk visc alpha ", alpha1) print("gas tau ",taugas1) print("depletion tau ",taudep1) print("gas gisk lifetime ",lifetime1) print("migration coefficient ", migratek1)

1. print (np.sum(ms1), np.sum(mall1))

len1=len(aas1)

1. print ("a_AU m_Me")

1. for n in range(0, len1):
2. print(round(aas1[n],3), round(ms1[n],4))

1. df1=pd.DataFrame([aas1, ms1])

df1=pd.DataFrame({'AU':aas1,'M':ms1,'R':rs1, "P":orbitperiods1, "ESI":esis1, "Ts/K":tsurfplanets1, "gravity": gees1, "vesc": vescs1, "Patm": patms1, "sol": sols1, "rotation": rots1, "locked": locked1,"minmolmass": minimol1, "magnetic": magnetic1,"rockp": rockp1, "icep": icep1, "gasp": gasp1, "type": planettypes1} )

df1=df1.sort_values(by=['AU'])

1. df1.columns=['a_AU', 'm_Me']

print() print("Generated planets: ") print("---------------------")

print(df1)

print("AU:") print(aas1) print("corresponding distance in solar system") print(correspondic_solar_distances1)

df1.to_csv(outfilename1, index=False)

1. quit(-1)

df1=df1.sort_values(by=["AU"])

aas2=df1['AU'].to_numpy() rs2=df1['R'].to_numpy() planettypes2=df1['type'].to_numpy()

generate_povray(name1, aas2, rs2, planettypes2, drawskala1)

plotting4(aas1, eccs1,ms1, 0, 40)

1. plt.show()

print(".")