2. 1. generate fictional solar system: masses, distances: a bit only solar system like, or not
4. code Python 3, Vapory POV-Ray rendering "sudo apt install povray && pip install vapory"
5. 1. testing code only
6. 1. 25.02.2024 0000.00011

import math import mpmath import numpy as np import random import sys import time

import matplotlib as mpl import matplotlib.pyplot as plt

1. python3 Pov-Ray extension 2024

from vapory import *

pi=math.pi degrad=pi/180

1. 1. si units

au=149597870700

G=6.6743e-11 kb=1.380649e-23 sb= 5.670374419e-8 planck=6.62607015e-34 avogadro=6.02214076e23 rydberg=10973731.6 mass_proton=1.67262192e-27

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

mearth=5.9722e24 rearth=6.3781e6 rooearth=5517 year1=31556926 megayear1=year1*1e6

mjupiter_me=317.828

mjupiter=1.89813e27

def experimental_gas_acc(mcorebegin1, sigmak1, time1, taudisk1):

mcore1=mcorebegin1*math.pow(sigmak1,1.5) ## exp maybe 2 or 1 ... 2.5 ## simplified and modified from ## https://arxiv.org/pdf/1904.10470.pdf # The Boundary Between Gas-rich and Gas-poor Planets # Draft version June 5, 2019 Eve J. Lee mgask_slow1=0.01*math.pow((sigmak1), 1/8)*math.pow((mcore1/10), 2)*math.pow((time1/0.1), 0.5) mgask_runaway1=math.exp((time1/taudisk1)*(math.pow((mcore1/10), 4)))

return(mcorebegin1, gas1)

def merge_planets_1(maxdeltak1, planetdistances0, planetmasses0, stones0, ices0, gases0, waters0): las1=planetdistances0.tolist() mas1=planetmasses0.tolist() sts1=stones0.tolist() ics1=ices0.tolist() gas1=gases0.tolist() wets1=waters0.tolist() las2=[] mas2=[] sts2=[] ics2=[] gas2=[] wets2=[] len1=len(las1) sweet1=0 m=1 for n in range(1, len1-2): a1=las1[m-1] a2=las1[m] m1=mas1[m-1] m2=mas1[m] s1=sts1[m-1] s2=sts1[m] i1=ics1[m-1] i2=ics1[m] g1=gas1[m-1] g2=gas1[m] w1=wets1[m-1] w2=wets1[m] a3=(a1+a2)/2 delta1=abs(a2-a1) mrel1=m1/(m1+m2) adda1=delta1*mrel1 #print(delta1, mrel1, adda1) #quit(-1) a3=a1+adda1 #print(sweet1, a3) #quit(-1) m3=(m1+m2) s3=(s1+s2) i3=(i1+i2) g3=(g1+g2) w3=(w1+w2)

maxdeltak2=(a2*maxdeltak1)

if(sweet1==0): #print("*")

if(delta1<maxdeltak2): #del las1[m] las2.append(a3) mas2.append(m3) sts2.append(s3) ics2.append(i3) gas2.append(g3) wets2.append(w3) sweet1=1 m=m+2 #del las1[m-2] len1=len1-1 print(las2) #quit(-1) else: las2.append(a1) mas2.append(m1) sts2.append(s1) ics2.append(i1) gas2.append(g1) wets2.append(w1) m=m+1 else: #print("¤¤") las2.append(a1) mas2.append(m1) sts2.append(s1) ics2.append(i1) gas2.append(g1) wets2.append(w1) m=m+1

las2.append(a2) mas2.append(m2) sts2.append(s2) ics2.append(i2) gas2.append(g2) wets2.append(w2)

las3=np.array(las2) mas3=np.array(mas2) sts3=np.array(sts2) ics3=np.array(ics2) gas3=np.array(gas2) wets3=np.array(wets2) return(las3,mas3,sts3, ics3, gas3, wets3)

def vapory_render_test_1(aas1, mas1, types1, tilts1, rings1, moons1):

   viewx1=200
   viewz1=800
   viewy1=0
   rcoeff1=0.05
   acoeff1=400
   viewangle1=60

   camera = Camera( 'location', [viewx1, viewy1, -viewz1], 'look_at', [viewx1, 0, 0] , "angle", viewangle1)
   light1 = LightSource( [-20000, 0, -20000], 'color', [1*1.5, 1*1.5, 1*1.5] )
   
   x1=[]
   rs1=[]
   colors1=[]
   len1=len(aas1)
   for n in range(0, len(aas1)-1):
       a1=aas1[n]
       r1=mas1[n]
       type1=planettypes1[n]
       ra1=r1*rcoeff1
       if(r1<5): ra1=ra1*300
       r1=ra1
   
       #print(a1, r1)
       x1.append(math.log10(a1)*acoeff1)
       rs1.append(r1)
       if(planettypes1[n]==0): colors1.append([1,.6,.5])
       if(planettypes1[n]==1): colors1.append([.5, .5, .9])
       if(planettypes1[n]==2): colors1.append([.8, 1, .2])

   objekts1=[]
   planets1=[]
   n=0
   objekts1.append(light1)
   
   #rings1[6]=1
   
   for n in range(0,(len(x1)+0)):
       type1=planettypes1[n]
       ringed1=rings1[n]
       tilt1=tilts1[n]
       rad1=rs1[n]
       xloc1=x1[n]
       #rin
       #sphere1 = Sphere( [x1[n], 0, 0], rs1[n], Texture( Pigment( 'color',colors1[n] ), Finish( 'phong', 0.05,'specular',0.01,'metallic', 0.0 ) )  )
       if(type1==0): ##rocky
           sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture(  Pigment( 'wrinkles', 'scale',[10,10,10], PigmentMap([0.00, 'rgb', [ 1,.6,.5]],[1.0 , 'rgb', [0.5,0.5,0.5]] ) ) , Finish( 'phong', 0.05,'specular',0.01,'metallic', 0.0 )) )           
         
       if(type1==4): ## terran
           sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture(  Pigment( 'function {f_wrinkles(x,y,z)*0.9+f_granite(x,y,z*0.1)}', 'scale',[30,30,30], PigmentMap( [0.00, 'rgb', [ 1,0.6,0.5]],[0.20, 'rgb', [ 0,1,0]],[0.5, 'rgb', [ 0,1,0]], [0.5 , 'rgb', [0,0,1]], [1.0 , 'rgb', [0,0,1]] ) ) , Finish( 'phong', 0.05,'specular',0.01,'metallic', 0.0 )), Texture(  Pigment( 'wrinkles', 'turbulence',2, 'scale',[20,20,20], PigmentMap([0.00, 'rgbt', [ 1,1,1,1]], [0.3, 'rgbt', [ 1,1,1,1]],[1.0 , 'rgbt', [1,1,1,0]] ) ) , Finish( 'phong', 0.1,'diffuse',0.6,'specular',0.01,'metallic', 0.0 )) )           

       if(type1==2): ## neptune icy subgiant or giant
           sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture(  Pigment( 'function  { y*2*f_wrinkles(5,y*3,5)*(1-f_wrinkles(x/5,y,z/5)) }', 'sine_wave turbulence', 0.1 , 'scale',[10,10,10], PigmentMap( [0.00, 'rgb', [ 0.3,0.3,0.8]],  [0.999, 'rgb', [ 0.4,0.4,0.9]],  [1.0 , 'rgb', [0.9,0.9,0.9]] ) ) , Finish('diffuse', 0.6  ) ), "rotate x*", tilt1)     
      
       if(type1==1): ## icy 
           sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture(  Pigment( 'function { f_wrinkles(x,y,z) }', 'turbulence', 0.05 , 'scale',[10,10,10], PigmentMap([0.00, 'rgb', [ 0.4,0.4,0.4]], [0.65 , 'rgb', [0.8,0.8, 0.8]] ,  [1.0 , 'rgb', [1,1,1]] ) ) , Finish('diffuse', 0.7  ) ) )     
       if(type1==10): ## saturn
           sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture(  Pigment( ' function { y*4*f_wrinkles(0,y*1,0) } sine_wave turbulence y*0.5 scale 20', 'turbulence', 0.02 , 'scale',[10,2,10], PigmentMap([0.00, 'rgb', [ 0.905882, 0.866667, 0.694118]], [1.0 , 'rgb', [0.870588, 0.752941, 0.611765]] ) ) , Finish('diffuse', 0.6 ) ) , "rotate x*", tilt1)  
       if(type1==3): ## gaseous giant
           sphere1=Sphere( [0, 0, 0], rad1, "rotate z*",tilt1, "translate x*", xloc1, Texture(  Pigment( 'function { ((sin((y+f_wrinkles(1,y,0))/2)+1)/2)*0.5+f_wrinkles(1,y/5,2)   }  sine_wave  turbulence y*0.05  scale 2  scale 0.2 warp {turbulence 0.2} scale 5 scale x*0', PigmentMap([0.00, 'rgb', [0.992157, 0.996078, 0.929412]], [1.0 , 'rgb', [ 0.984314, 0.721569, 0.52549]] ,  [1.0 , 'rgb', [ 0.368627, 0.333333, 0.313725 ]]) ) , Finish( 'diffuse', 0.6 ) ), "rotate x*", tilt1 )  
           
       if(type1==5): ## venusian
           sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture(  Pigment( 'wrinkles turbulence 0.4', 'scale',[10,10,10], PigmentMap([0.00, 'rgb', [ 0.423529, 0.278431, 0.2]],[1.0 , 'rgb', [0.917647, 0.917647, 0.917647]] ) ) , Finish( 'phong', 0.05,'specular',0.01,'metallic', 0.0 )) )    
       if(type1==6): ## mercurial
           sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture(  Pigment( 'granite', 'scale',[10,10,10], PigmentMap([0.00, 'rgb', [ 1,1,1]],[1.0 , 'rgb', [0.5, 0.5, 0.5]] ) ) , Finish( 'phong', 0.05,'specular',0.01,'metallic', 0.0 )) )    
       if(type1==7): ## asteroids
           sphere1=Sphere( [x1[n], 0, 0], 30, Texture(  Pigment( 'granite', 'turbulence',0, 'scale',[50,50,15], PigmentMap([0.00, 'rgbt', [ 1,1,1,0]],[0.1 , 'rgbt', [1,1,1,0]],  [0.1 , 'rgbt', [1,1,1,1]], [1.0 , 'rgbt', [1,1,1,1]] ) ) )  )    
       if(type1==8): ## ocean
           sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture(  Pigment( 'granite', 'turbulence',1, 'scale',[20,20,20], PigmentMap([0.00, 'rgbt', [ 0,0,1,0]], [1.0 , 'rgbt', [1,1,1,0]] ) ) )  )
       if(type1==9): ## cold stony
           sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture(  Pigment( 'wrinkles', 'turbulence',0.0, 'scale',[10,10,10], PigmentMap([0.00, 'rgbt', [ 1,.6,.5,0]], [1.0 , 'rgbt', [1,1,1,0]] ) ) )  )
       if(ringed1==1): ## ringed1
           #ring1=Torus([0,1],"scale <2,0.1,2>*",rad1, "rotate x*45 translate x*",xloc1)

1. ring1=Sphere([0,0,0], rad1*2.5 ,"scale y/200 rotate x*90 rotate z*0","translate x*", xloc1 , Texture( Pigment( 'function {f_wrinkles(sqrt(x*x+y*y+z*z),sqrt(x*x+y*y+z*z) ,sqrt(x*x+y*y+z*z) ) }', ' rotate x*90 turbulence',0.0, 'scale',[1,1,1], PigmentMap([0.00, 'rgbt', [ 2,1,1,0]], [1.0 , 'rgbt', [1,1,1,0]] ) ) ) )
2. ring1=Sphere([0,0,0], rad1*2.5 ,"scale y/200 rotate x*90 rotate z*0","translate x*", xloc1 , Texture( Pigment( 'function {f_onion(x,y,z )*f_wrinkles(sqrt(x*x+y*y+z*z),1,2) }', ' rotate x*00 turbulence',0.0, 'scale',[20,20,20], "translate x*",xloc1, PigmentMap([0.00, 'rgbt', [ 1,1,1,0]], [1.0 , 'rgbt', [1,1,1,1]] ) ) ), "rotate x*",tilt1 )
3. ring1=Torus(rad1*2 ,rad1*0.75,"scale y/200 rotate x*90 rotate z*0","translate x*", xloc1 , Texture( Pigment( 'function {f_onion(f_wrinkles1(x,y,z),y,z )*f_wrinkles(sqrt(x*x+y*y+z*z),1,2) }', ' rotate x*00 turbulence',0.0, 'scale',[20,20,20], "translate x*",xloc1, PigmentMap([0.00, 'rgbt', [ 1,1,1,0]], [1.0 , 'rgbt', [1,1,1,1]] ) ) ), "rotate x*",tilt1 )

           ring1=Torus(rad1*1.5 ,rad1*0.33,"scale y/200 rotate x*0 rotate z*0","translate x*", xloc1 , Texture(  Pigment( 'function { f_granite(f_onion(x,y,z),2,3) }', 'sine_wave rotate x*00 turbulence',0.0, 'scale',[20,20,20], "translate x*",xloc1, PigmentMap([0.00, 'rgbt', [ 1,1,1,1]], [1.0 , 'rgbt', [1,1,1,0]] ) ) ), "rotate x*",tilt1  )

           objekts1.append(ring1)

1. if(type1==9): ## cold stony
2. sphere1=Sphere( [x1[n], 0, 0], rs1[n], Texture( Pigment( 'wrinkles', 'turbulence',0.0, 'scale',[10,10,10], PigmentMap([0.00, 'rgbt', [ 1,.6,.5,0]], [1.0 , 'rgbt', [1,1,1,0]] ) ) ) )

#if(ringed==1): # ring1=Torus( [0,1], rs1[n] "translate x*",x1, Texture( Pigment( Finish( 'diffuse',0.8,'phong', 0.05,'specular',0.01,'metallic', 0.0 )) ) # objekts1.append(ring1)

       planets1.append(sphere1)

       objekts1.append(planets1[n])
   
   scene1 = Scene( camera, objects= objekts1,   included=['colors.inc', 'textures.inc', 'stones.inc', 'skies.inc', 'functions.inc'],
             global_settings=['assumed_gamma 2.2'])
   scene1.render("lanets1", width=1600, height=200 , antialiasing=0.3, quality=9)    

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 = pow(((solarconstant*(1-albedo)) / (4*stefan*emission_ir)), 1/4)
   return(tee1)

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=math.sqrt(((rau2*rau2*rau2)/massstar)) Q1=100 ## earth today 13, can be 100

rt=0.027*math.pow( ((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 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 rotation_period_mr(mass, radius): # https://arxiv.org/pdf/2301.13836.pdf rotdays=0.632*math.pow(mass, 1/2)*radius # rotdays= 0.632*math.pow(mass, 1/6)*math.pow(rooearths, 1/3) return(rotdays)

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 new_gas_add1(a1,snowline1, isomass1):

   q1=a1/snowline1 ## or forst giant loc giantloc1
   gasadd0=310*math.pow(q1, -2) ## -2-2 desch
   gasadd1=gasadd1-isomass1     
   return(gasadd1)  

def merge_orbits_1(maxdeltak1, planetdistances1):

   len1=len(planetdistances1)
   ad2=[]
   md2=[]

   a1=planetdistances1[0]
   #m1=planetmasses1[0]

   for n in range(1, len1):
       a2=planetdistances1[n]
       #m2= planetmasses1[n]
       delta1=abs(a2-a1)
       maxdelta1=(1*maxdeltak1)*a2
       print(delta1, maxdelta1)
       if(delta1<maxdelta1):
           print("abby")
           a3=(a2+a1)/2
           ad2.append(a3)
           #m3=m1+m2
           #a1=a3
           #m1=m3
           n=n+2
       else:
           #m1=m2
           a3=a1
           a1=a2            
           ad2.append(a3)
      
       #md2.append(m1)

   #planetdistances2=planetdistances1
   planetdistances2=np.array(ad2)
   #planetmasses2=np.array(md2)
   print(planetdistances1)
   print(planetdistances2)
   #quit(-1)
   return(planetdistances2)

def current_milli_time():

   return round(time.time() * 1000)

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 water_content_dryplanets(aus1, mas1, snowline1): ## giant planet filters water, if just "above snowline wetness0=(aus1-snowline1)/snowline1 #wetness1=np.power(wetness0, 4) #wetness1=np.power(np.exp(wetness0),2) #wetness1=np.power(wetness0, 4) wetness2=0.1/np.sqrt(2*math.pi)*np.exp(-(wetness0*wetness0)/2) wetness1=np.power(wetness2,4)*1000*mas1 wetness1=np.where(wetness1>0.0015, 0.1, wetness1) wetness1=np.where(wetness1>0.75, 0.75, wetness1) return(wetness1)

def water_content_wetplanets(aus1, mas1, snowline1): ## non water filtering of giant planet wetness0=(np.copy(aus1))/snowline1 wetness1=(wetness0) wetness1=np.where(wetness1>0.75, 0.75, wetness1) return(wetness1)

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)
   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 migrate_random_inward(logamount1, numplanets1, aas1, ms1):

   migrate1=math.pow(10, -logamount1)
   #print(migrate1)
   migrate2=(0.5+np.random.rand()*0.5)*migrate1*10   
   ms2=ms1
   #aas2=aas1*1/math.exp(np.random.rand())
   #aas2=aas1*math.pow(np.random.rand(),3)
   aas2=aas1*migrate2
   return(aas2, ms2)

def random_gas_mass_1(a1, giantloc1, addgas1):

   giantcoeff1=2*giantloc1
   q1=abs(np.random.normal(loc=0, scale=1, size=len(a1)))
   p2=np.exp(-(a1-giantloc1)/giantcoeff1)
   #    p2=math.pow((a1-giantloc1)/(giantloc1*3), -1)
   q2=addgas1*q1*p2
   return(q2)    

def create_planet_distances_2(planetnumber1):

   #ps1=np.random.normal(loc=2, scale=1, size=planetnumber1).astype(int)
   #psd1=np.random.normal(loc=1, scale=1, size=planetnumber1).astype(int)+1
   ps1=np.random.normal(4,size=planetnumber1).astype(int)
   psd1=np.random.normal(4, size=planetnumber1).astype(int)
   ps2=ps1+psd1            
   print(ps1)
   print(ps2)
   pr1=ps2/ps1
   #print(pr1)
   pr2=pr1*0
   len1=len(pr1)
   pera1=1/4
   per1=pera1
   for n in range(0, len1):
       pr2[n]=per1
       per1=per1*pr1[n]

   #print(pr2)
   asa1=np.power(pr2, 2/3)
   #print(asa1)
   asa2=asa1+asa1*np.random.normal(loc=0, scale=0.05, size=planetnumber1)
   print(asa2)
   #quit(-1)            
   return(asa2)

def create_planet_distances_1(planetnumber1):

   distance_base1=0.3+np.random.normal(loc=0, scale=0.2)
   distance_coeff1=0.3+np.random.normal(loc=0, scale=0.1)
   distance_exp1=1.8+np.random.normal(loc=0, scale=0.4)

   distans1=np.linspace(1, planetnumber1, planetnumber1+1)
   disnansmover1=np.random.normal(loc=1, scale=0.3)
   distans1=np.power(distance_base1+(distans1)*distance_coeff1, distance_exp1)
   #print(distans1)
   #quit(-1)
   return(distans1)

def create_planet_system(seed1, starmass1, metallicity1, gassy1, disktau1):

   #seed1=123
   discecc1=0.1 ## not used
   starlum1=math.pow(starmass1, 3.5) ## meybe nok, because pre-ms lum is diffenent
   random.seed(seed1)
   np.random.seed(seed1)
   planetnumber1=10
   snowline1=5.0*starmass1 ## maybe nok different distance in pre-ms
   snowline2=1*math.sqrt(starlum1) ## main sequence inner snowy line
   kuiperlimit1=45*starmass1
   icecoeff1=4
   gaslimit1=8
   maxisomass1=25 #*math.pow(starmass1, 1/6)
   giant_effectlimit=30 ## planet can create inner asteroid belt
   massk1=1*np.random.normal(loc=1.2, scale=0.5) #*math.pow(starmass1, 1/6)
   oligarch_collisions1=1
   gasamount_coeff1=500*math.exp(disktau1/2.8)*gassy1*math.pow(starmass1*metallicity1, 1/6)
   radiusc1=25*starmass1 ## disk edge
   radiusc2=5*starmass1 ## disk edge sharpness
   #gasadd_distance_exponent1=-3/1
   gasadd_distance_exponent1=-3/1
   planetdistances1=create_planet_distances_2(planetnumber1+1)
   #planetdistances1=create_planet_distances_1(planetnumber1+1) 
   if(planetdistances1.any()==np.nan):
           planetdistances1=create_planet_distances_1(planetnumber1+1)        
   if(planetdistances1.any()==np.inf):
           planetdistances1=create_planet_distances_1(planetnumber1+1) 

   planetdistances1=planetdistances1*math.sqrt(starlum1)
   deltak1=2
   
   
   #planetdistances1=merge_orbits_1(deltak1, planetdistances1)
   #planetnumber1=len(planetdistances1)-1

   print("Planets # ", planetnumber1)
   
   snowy1=np.copy(planetdistances1)
   snowy1=np.where(planetdistances1>snowline1,icecoeff1,1)
   icy1=np.where(planetdistances1>snowline1,1,0)
   planetmasses1=starmass1*metallicity1*massk1*0.1*np.power(planetdistances1,3)*snowy1*np.power(planetdistances1,-3/2)*np.power(planetdistances1/10,-1)
   #planetmasses1=planetmasses1*np.random.randint(4)
   planetmasses1=planetmasses1*np.random.random(planetnumber1+1)*oligarch_collisions1*(0.75+np.random.rand(planetnumber1+1)*0.5) ## 4 origarch collisions
   planetmasses1=np.where( planetmasses1>maxisomass1, maxisomass1, planetmasses1)        
   ices1=planetmasses1*icy1*(icecoeff1-1)/icecoeff1
   stones1=planetmasses1/icecoeff1
   snowy1=(np.copy(planetdistances1)/snowline1)
   snowy1=np.where(snowy1>1,1,0)
   gasmassy1=(np.copy(planetmasses1)/gaslimit1).astype(int)
   gasadder1=gasmassy1*snowy1
   gasamount1=random_gas_mass_1(planetdistances1, snowline1, gasamount_coeff1)*gasadder1
   densitycoeff1=1-np.exp( (planetdistances1-radiusc1)/radiusc2 )/1000
   densitycoeff1=np.where( densitycoeff1<0, 0, densitycoeff1)
   ## also add rocks and ice with gas!
   iceadd1=gasadder1*0.75*0.1
   stoneadd1=gasadder1*0.25*0.1
   ices1=ices1+iceadd1
   stones1=stones1+stoneadd1
   planetmasses3=planetmasses1+gasamount1+iceadd1+stoneadd1
   giants1=np.argwhere(planetmasses3>giant_effectlimit)
   ## giants at snow line reduce mass of inner disk!

   #wetness1=np.power((planetdistances1/snowline1),4)/5*planetmasses3
   ## default: inner wetness1
   wetness1=water_content_wetplanets(planetdistances1, planetmasses3, snowline2)

   #wetness1[1]=0  

   #print(giants1)
   if len(giants1>0):
       #wetness1=np.power((planetdistances1/snowline1),4)*planetmasses3
       wetness1=water_content_dryplanets(planetdistances1, planetmasses3, snowline1)
       #wetness1[0]=0  
       giant=int(giants1[0])
       #print(giant)
 
       # giants reduces mass of inner planets

       if(giant>0):
           giant_distance1=planetdistances1[giant]
           giant_mass1=planetmasses1[giant]
           #planetmasses3[giant-1]=0.001   
           for n in range(0, giant):
               diff1=planetdistances1[n]/giant_distance1
               au1=planetdistances1[n]
               if (diff1>0.4):planetmasses3[n]=0.01 ## asteroids
               if (diff1>0.25):planetmasses3[n]=planetmasses3[n]/10 ## mars-like
               #if (diff1<0.1):planetmasses3[n]=planetmasses3[n]/20  ## mercury-like
               if (au1<0.5):planetmasses3[n]=planetmasses3[n]/20  ## mercury-like
               
               #print(diff1)

   wetness1=np.where(wetness1>0.75, 0.75,wetness1)
   #print(wetness1)

   ## Kuiper belt!
   for n in range(0, planetnumber1):
       if(planetdistances1[n]>kuiperlimit1):
           planetmasses3[n]=1e-3
           ices1[n]=planetmasses3[n]*3/4
           stones1[n]=planetmasses3[n]*1/4            
           gasamount1[n]=0       
           wetness1[n]=0.75

   return(planetdistances1, planetmasses3, stones1, ices1, gasamount1, wetness1)

1. 1. 1. 1. 2. main program

starmass1=1 metallicity1=1.0 ##solar coeff abs, not log gasamount1=1.0 ## disktau1=2.0 ## diskage1=1.0 ## not used

starlum1=math.pow(starmass1, 3.5)

1. print(rad1)

1. quit(-1)

1. seed1= 73255405
2. seed1= 1564709
3. seed1=321

1. 1. random seed:
2. seed1=None
3. seed1 = int(random.randrange(sys.maxsize)/1e12)

seed1=current_milli_time()%100000000

1. seed1=1

1. seed1=321

1. seed1= 89881932

1. seed1= 1564709

print("Seed is: ", seed1)

1. rnd1 = random.Random()
2. seed1 = rnd1.get_seed()

1. print("Seed is ", seed1)

1. create_planet_distances_2(10)

1. quit(-1)

for mama in range(0,100): seed1=current_milli_time()%100000000 #seed1= 89881932 seed1= 1564709 #seed1= 92365108 #seed1= 92710846 print("Seed is ", seed1) planetdistances1, planetmasses1, stones1, ices1, gases1, wetness1=create_planet_system(seed1, starmass1, metallicity1, gasamount1, disktau1) maxdeltak1=0.33 ## Testing ... planetdistances1, planetmasses1, stones1, ices1, gases1, wetness1= merge_planets_1(maxdeltak1, planetdistances1, planetmasses1, stones1, ices1, gases1, wetness1) planetnumber1=len(planetdistances1) #logamount1=2 #planetdistances1, planetmasses1 =migrate_random_inward(logamount1, planetnumber1, planetdistances1, planetmasses1) planetmosses1=stones1 + ices1 + gases1 stonepercent1=(100*stones1/planetmosses1) icepercent1=(100*ices1/planetmosses1) gaspercent1=(100*gases1/planetmosses1)

planetradiuses1=np.power(planetmasses1,1/3) for n in range(0, len( planetmasses1)): planetradiuses1[n]=radius_s1(planetmasses1[n], stonepercent1[n], icepercent1[n], gaspercent1[n])

planettemps1=planet_temperatures_1(starmass1, starlum1, planetdistances1) planetesis1=calc_esis(planetmasses1, planetradiuses1,planettemps1)

planettypes1=np.copy(planetdistances1)*0

planettilts1=abs(np.random.normal(loc=0, scale=20, size=len(planetdistances1))) planetmoons1=planettilts1*0 planetrings1=(np.copy(planetmasses1)/20).astype(int) planetrings1=np.where(planetrings1>0,1,0)*np.random.randint(2, size=len(planetrings1))

#print(planetrings1) #quit(-1)

for n in range(0, len( planetmasses1)): #default stone, rock planettypes1[n]=0 ## rocky if(icepercent1[n]>1): planettypes1[n]=2 ## icy thick atm if(icepercent1[n]>1): if(planetmasses1[n]<4): ## icy thin atm planettypes1[n]=1 if(gaspercent1[n]>1): ## gaseous giant planettypes1[n]=3

if(gaspercent1[n]>1): ## gaseous "saturn" if(planetmasses1[n]<150): if(planetmasses1[n]>30): planettypes1[n]=10 if(stonepercent1[n]>99.9): ## almost pure rocky if(planetesis1[n]>0.95): planettypes1[n]=4 ## terran ESI deviation 5%

if(stonepercent1[n]>99.9): ## almost pure rocky if(planettemps1[n]>340): if(planetmasses1[n]>0.3): planettypes1[n]=5 ## venusian if(planetmasses1[n]<=0.3): planettypes1[n]=6 ## mercurial if(planettemps1[n]<=240): if(stonepercent1[n]>99.9): ## almost pure rocky planettypes1[n]=9 ## icy stony if(planettemps1[n]>=273): if(planettemps1[n]<=373): if(icepercent1[n]>0.1): ## almost pure rocky if(gaspercent1[n]<1): ## almost pure rocky planettypes1[n]=8 ## ocean stony print("OCEAN") if(planetmasses1[n]<=0.001): planettypes1[n]=7 ## asteroids if (planetesis1.any()>0.95): if(planetmasses1.any()>300): print("Tries ", mama) break

planeteccs1=abs(np.random.normal(loc=0, scale=0.05, size=len(planetdistances1))) planetorbincs1=abs(np.random.normal(loc=0, scale=1, size=len(planetdistances1))) planetperihels1=np.random.random(size=len(planetdistances1))*360

planetorbitperiods1=np.power(planetdistances1, 3/2)

1. print(planetorbitperiods1)

1. quit(-1)

planetsols1=np.power(planetdistances1, -2)*math.sqrt(starlum1) planetdensities1=planetmasses1/np.power(planetradiuses1,3) planetvescs1=np.sqrt(planetmasses1/planetradiuses1) planetgravities1=planetmasses1/np.power(planetradiuses1,2)

planetatmpressures1=atmosphere_pressure_estimation(planetmasses1, planetradiuses1,planettemps1)

1. 1. simple approx

planetrotperiods1=24*np.power(planetmasses1, -1/8)+1*np.random.normal(loc=0, scale=1, size=len(planetdistances1))

1. 1. oligarch collision changes tilt, maybe rot period

planetmagnetic1=magnetic_potential_estimation(planetmasses1, planetradiuses1, planetrotperiods1)

1. 1. stony object, big mass slower rot

for n in range(0, len( planetmasses1)):

   if(planetdensities1[n]>0.8):
       planetrotperiods1[n]=24*np.power(planetmasses1[n], -1/8)

1. 1. simple resonance locked approx

for n in range(0, len( planetdistances1)):

   if(planetdistances1[n]<0.25): planetrotperiods1[n]=24*365*planetorbitperiods1[n]
   if(planetdistances1[n]<0.5): planetrotperiods1[n]=24*365*planetorbitperiods1[n]*2/3
   if(planetdistances1[n]<0.75): planetrotperiods1[n]=24*365*planetorbitperiods1[n]*3/4

1. initial albedos

planetalbedos1=planetdistances1*0+0.2 for n in range(0, len( planetmasses1)):

   planetalbedos1[n]=0.30 #default
   if(planettypes1[n]==4): ## terran
       planetalbedos1[n]=0.31
   if(planettypes1[n]==0): ## generic, desert
       planetalbedos1[n]=0.2
   if(planettypes1[n]==1): ## icy
       planetalbedos1[n]=0.76
   if(planettypes1[n]==2): ## neptune
       planetalbedos1[n]=0.32               
   if(planettypes1[n]==2): ## jupiter
       planetalbedos1[n]=0.51   
   if(planettypes1[n]==10): ## saturn
       planetalbedos1[n]=0.51   
   if(planettypes1[n]==7): ## asteroids
       planetalbedos1[n]=0.1   
   if(planettypes1[n]==6): ## moon, mercurial
       planetalbedos1[n]=0.09   

1. print(planetdistances1)
2. print(planettemps1)
3. print(planetesis1)

1. print(magnetic1)

1. quit(-1)

print ("planet# distance mass radius stone% ice% gas% Temp./K ESI ") for n in range(0, planetnumber1, 1):

   print(str(n+1).rjust(3), str(round(planetdistances1[n],2)).rjust(10) , str(round(planetmasses1[n],3)).rjust(10) ,  str(round(planetradiuses1[n],3)).rjust(10) ,str(round(stonepercent1[n],1)).rjust(10),  str(round(icepercent1[n],3)).rjust(10), str(round(gaspercent1[n],3) ).rjust(10) ,str(round(planettemps1[n],3)).rjust(10), str(round(planetesis1[n],3)).rjust(10) )

print ("planet# density gravity vesc P_atm magnetic_pot rotperiod_h ecc axis_tilt") for n in range(0, planetnumber1, 1):

   print(str(n+1).rjust(3), str(round(planetdensities1[n],2)).rjust(10), str(round(planetgravities1[n],2)).rjust(10),str(round(planetvescs1[n],2)).rjust(10),str(round(planetatmpressures1[n],2)).rjust(10), str(round(planetmagnetic1[n],2)).rjust(10),  str(round(planetrotperiods1[n],2)).rjust(10),   str(round(planeteccs1[n],2)).rjust(10),  str(round(planettilts1[n],2)).rjust(10),  )

print ("planet# density_gcm2 albedo") for n in range(0, planetnumber1, 1):

   print(str(n+1).rjust(3), str(round(planetdensities1[n]*5.6,2)).rjust(10) ,  str(round(planetalbedos1[n]*1,2)).rjust(10)  )

stdout0 = sys.stdout fout1 = open('planets.txt', 'w') sys.stdout = fout1

print ("planet# distance mass radius stone% ice% gas% Temp./K ESI ") for n in range(0, planetnumber1, 1):

   print(str(n+1).rjust(3), str(round(planetdistances1[n],2)).rjust(10) , str(round(planetmasses1[n],3)).rjust(10) ,  str(round(planetradiuses1[n],3)).rjust(10) ,str(round(stonepercent1[n],1)).rjust(10),  str(round(icepercent1[n],3)).rjust(10), str(round(gaspercent1[n],3) ).rjust(10) ,str(round(planettemps1[n],3)).rjust(10), str(round(planetesis1[n],3)).rjust(10) )

print ("planet# density gravity vesc P_atm magnetic_pot rotperiod_h ecc axis_tilt") for n in range(0, planetnumber1, 1):

   print(str(n+1).rjust(3), str(round(planetdensities1[n],2)).rjust(10), str(round(planetgravities1[n],2)).rjust(10),str(round(planetvescs1[n],2)).rjust(10),str(round(planetatmpressures1[n],2)).rjust(10), str(round(planetmagnetic1[n],2)).rjust(10),  str(round(planetrotperiods1[n],2)).rjust(10),   str(round(planeteccs1[n],2)).rjust(10),  str(round(planettilts1[n],2)).rjust(10),  )

print ("planet# density_gcm2 albedo") for n in range(0, planetnumber1, 1):

   print(str(n+1).rjust(3), str(round(planetdensities1[n]*5.6,2)).rjust(10) ,  str(round(planetalbedos1[n]*1,2)).rjust(10)  )

sys.stdout = stdout0 fout1.close()

1. for x in range(1, 11):
2. print('{0:2d} {1:3d} {2:4d}'.format(x, x*x, x*x*x))

print (planetdistances1) print (planetmasses1)

print("seed1=", seed1)

vapory_render_test_1(planetdistances1, planetmasses1, planettypes1, planettilts1, planetrings1, planetmoons1)

1. quit(-1)

1. pgiants1<-np.copy()

1. plt.scatter(planetdistances1, planetmasses3*0+5, s=planetmasses3*30)

plt.gca().set_yticks([]) plt.yticks([])

1. ax.set_yticks([])

for m in range(0,planetnumber1): sx1=planetradiuses1[m]*100 plt.scatter(planetdistances1[m], 5, s=sx1, facecolor="#ff7030", color="#7f3010") if(planettypes1[m]==1): plt.scatter(planetdistances1[m], 5, s=sx1, facecolor="#dfdfff", color="#3f3fff") if(planettypes1[m]==2): plt.scatter(planetdistances1[m], 5, s=sx1, facecolor="#afafff", color="#3f3fff") if(planettypes1[m]==3): plt.scatter(planetdistances1[m], 5, s=sx1, facecolor="#a0ffa0", color="#007f00") if(planettypes1[m]==10): plt.scatter(planetdistances1[m], 5, s=sx1, facecolor="#c0ffc0", color="#00af00") if(planettypes1[m]==4): plt.scatter(planetdistances1[m], 5, s=sx1, facecolor="#a0a0ff", color="#a07f10") plt.scatter(planetdistances1[m], 5, s=sx1*0.3, facecolor="#a7ffa7", color="#00ff00")

plt.xlim(0,10)

1. 1. time when is mass recahed mnp.power(a, 4) or np.power(a, 4.2) a(3.6)
   2. 2x mmsn math.pow(a,6.3)

1. plt.plot(aas1, sigmas1)
2. plt.xlim(0.5, 30)
3. plt.ylim(0, 30)

1. plt.plot(aas1, isomass1)
2. plt.plot(aas1, growthrates1)

1. plt.plot(aas1, growthrates1*isomass1)

plt.show()