File talk:Planetary interior structure.png
Jump to navigation
Jump to search
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import time
import winsound
import math
global ED
ED=5.515 # Earth ElementDensity in g/cm3
JM=317.828 # Jupiter mass in Earth masses
JR=11.21 # Jupiter Radius in Earth radii
Elements=5
Nl=10 # Number of degeneracy levels
Mal=5 # Number of mass addition levels. This same variable is later resused to indicate which level it is on
CoreTemperature= np.zeros(Nl+2, dtype=np.float64)
BindingEnergy= np.zeros(Elements*Nl+2, dtype=np.float64)
Count= np.zeros(Elements*Nl+2, dtype=np.int)
Delta= np.zeros(Elements*Nl+2, dtype=np.float64)
DP= np.zeros(Elements*Nl+2, dtype=np.float64)
Element= np.zeros(Elements*Nl+2, dtype=np.int)
Finished= np.zeros(Elements*Nl+2, dtype=np.int)
GhostMass= np.zeros(Elements*Nl+2, dtype=np.float64)
Level= np.zeros(Elements*Nl+2, dtype=np.int)
Density= np.zeros(Elements*Nl+2, dtype=np.float64)
Volume= np.zeros(Elements*Nl+2, dtype=np.float64)
Mass= np.zeros(Elements*Nl+2, dtype=np.float64)
OriginalLayer= np.zeros(Elements*Nl+2, dtype=np.int)
PostCompressionLayer=np.zeros(Elements*Nl+2, dtype=np.int)
Pressure= np.zeros(Elements*Nl+2, dtype=np.float64)
Radius= np.zeros(Elements*Nl+2, dtype=np.float64)
TotalBindingEnergy= np.zeros(Elements*Nl+2, dtype=np.float64)
TotalMass= np.zeros(Elements*Nl+2, dtype=np.float64)
TotalPressure= np.zeros(Elements*Nl+2, dtype=np.float64)
TotalVolume= np.zeros(Elements*Nl+2, dtype=np.float64)
CuriePoint= np.zeros((Elements+2, Nl+2), dtype=np.float64)
BindingEnery= np.zeros((Elements+2, Nl+2), dtype=np.float64)
Degeneracy_Pressure= np.zeros((Elements+2, Nl+2), dtype=np.float64)
ElementDensity= np.zeros((Elements+2, Nl+2), dtype=np.float64)
Layer= np.zeros((Elements+2, Nl+2), dtype=np.int)
Delta2= np.zeros((Mal+1, Elements+2, Nl+2), dtype=np.float64)
Ratio= np.zeros((Elements*Nl+2, Elements*Nl+2), dtype=np.float64)
Hydrogen=1
Methane=2
Diamond=2
Water=3
Rock=4
Silicon=4
Iron=5
P=1.1 # increase in density due to pressure
HD=0.125
S=(2**3)
IS=1/S
ElementDensity[Hydrogen,1]=P*HD/ED # Pure helium atmosphere
ElementDensity[Hydrogen,2]=P*S*S*0.079/ED #
ElementDensity[Hydrogen,3]= S* ElementDensity[Hydrogen,2] #
ElementDensity[Hydrogen,4]= S* ElementDensity[Hydrogen,3] #
ElementDensity[Hydrogen,5]= S* ElementDensity[Hydrogen,4] #
ElementDensity[Methane,1]= P*0.422/ED #
ElementDensity[Diamond,2]= P*3.5/ED #
ElementDensity[Diamond,3]= P*HD*3*(3**3)/ED # ground shell
ElementDensity[Diamond,4]= P*S*S*ElementDensity[Diamond,3] #
ElementDensity[Diamond,5]= S* ElementDensity[Diamond,4] #
ElementDensity[Water,1]= P*1/ED #
ElementDensity[Water,2]= P*HD*IS*4*(4**3)/ED #
ElementDensity[Water,3]= S* ElementDensity[Water,2] # ground shell
ElementDensity[Water,4]= S*S*ElementDensity[Water,3] #
ElementDensity[Water,5]= S* ElementDensity[Water,4] #
ElementDensity[Rock,1]= P*3.346/ED #
ElementDensity[Rock,2]= P*HD*IS*6.5*(6.5**3)/ED #
ElementDensity[Rock,3]= S* ElementDensity[Rock,2] # ground shell
ElementDensity[Rock,4]= S*S*ElementDensity[Rock,3] #
ElementDensity[Rock,5]= S* ElementDensity[Rock,4] #
ElementDensity[Iron,1]= P*6.3/ED #
ElementDensity[Iron,2]= S* ElementDensity[Iron,1] #
ElementDensity[Iron,3]= S* ElementDensity[Iron,2] #
ElementDensity[Iron,4]= P*HD*14*(13**3)/ED # ground shell
ElementDensity[Iron,5]= S*S*ElementDensity[Iron,4] #
for x in range(1, Elements+1, 1):
for y in range(6, Nl+1, 1):
ElementDensity[x,y]= S*ElementDensity[x,y-1]
x=4*3.14159265*0.333333
EarthVolume=x*(1216+2270+2885)**3
CoreVolume=x*(1216+2270)**3
InnerCoreVolume=x*(1216)**3
EM = 5.9722*10**24 # Kg
ED=10**-12*EM/EarthVolume
CoreVolume=CoreVolume/EarthVolume
InnerCoreVolume=InnerCoreVolume/EarthVolume
EarthVolume=1
EM=1
MantleVolume=EarthVolume-CoreVolume
OuterCoreVolume=CoreVolume-InnerCoreVolume
InnerCoreMass=(EM - (MantleVolume*ElementDensity[Rock,1])-(OuterCoreVolume*ElementDensity[Iron,1]) )
InnerCoreDensity=InnerCoreMass/(InnerCoreVolume)
InnerInnerCoreVolume=InnerCoreVolume*(InnerCoreDensity-ElementDensity[Iron,2])/(ElementDensity[Iron,3]-ElementDensity[Iron,2])
VolumeOfIron = OuterCoreVolume + (InnerCoreVolume-InnerInnerCoreVolume)*ElementDensity[Iron,2]/ElementDensity[Iron,1] + InnerInnerCoreVolume*ElementDensity[Iron,3]/ElementDensity[Iron,1]
#raise SystemExit()
InitialRockVolume=0
InitialIronVolume=0
MantleVolume=(MantleVolume-InitialRockVolume)
VolumeOfIron=(VolumeOfIron-InitialIronVolume)
Delta2[1, 0, 0]=0
Delta2[1, Hydrogen,1]=0
Delta2[1, Methane,1]=0 # amount of methane
Delta2[1, Water,1]=0 # amount of water
Delta2[1, Rock,1]=MantleVolume # amount of rock to add each time
Delta2[1, Iron,1]=VolumeOfIron # amount of iron to add each time
Delta2[2, 0, 0]=3.6
Delta2[2, Hydrogen,1]=0
Delta2[2, Methane,1]=0 # amount of methane
Delta2[2, Water,1]=50 # amount of water
Delta2[2, Rock,1]=MantleVolume # amount of rock to add each time
Delta2[2, Iron,1]=VolumeOfIron # amount of iron to add each time
Delta2[3, 0, 0]=20
Delta2[3, Hydrogen,1]=0
Delta2[3, Methane,1]=75 # amount of methane
Delta2[3, Water,1]=50 # amount of water
Delta2[3, Rock,1]=MantleVolume # amount of rock to add each time
Delta2[3, Iron,1]=VolumeOfIron # amount of iron to add each time
Delta2[4, 0, 0]=40
Delta2[4, Hydrogen,1]=0
Delta2[4, Methane,1]=75 # amount of methane
Delta2[4, Water,1]=50 # amount of water
Delta2[4, Rock,1]=MantleVolume # amount of rock to add each time
Delta2[4, Iron,1]=VolumeOfIron # amount of iron to add each time
Delta2[5, 0, 0]=0.8*JM
Delta2[5, Hydrogen,1]=11715
Delta2[5, Methane,1]=75 # amount of methane
Delta2[5, Water,1]=50 # amount of water
Delta2[5, Rock,1]=MantleVolume # amount of rock to add each time
Delta2[5, Iron,1]=VolumeOfIron # amount of iron to add each time
CoreTemperature[1]=1.0
CoreTemperature[2]=5.0
CoreTemperature[3]=5.0
for x in range(4, Nl+1, 1):
CoreTemperature[x]=7.5
CuriePoint[Iron,1] = 1.1 # 1000 kelvins
CuriePoint[Water,3] = 6.0 # 6000 kelvins
HeliumRadius=1.8506
CarbonRadius=1.2647/HeliumRadius
OxygenRadius=2.025/HeliumRadius
IronRadius =1.2448/HeliumRadius
HeliumRadius=1
HDP = (4.9*10**6)/(3.45*10**6) # helium degeneracy pressure
SDP = 4 * HDP # Standard degeneracy pressure
MF=1/32 # Metal factor
MMF=0.25 # Magnetic metal factor. Only affects iron at low temperature
P=6
S=2**P
Degeneracy_Pressure[Hydrogen,1] = HDP # Ground shell
Degeneracy_Pressure[Hydrogen,2] = S * S * MF * Degeneracy_Pressure[Hydrogen,1]
Degeneracy_Pressure[Hydrogen,3] = S * Degeneracy_Pressure[Hydrogen,2]
Degeneracy_Pressure[Hydrogen,4] = S * Degeneracy_Pressure[Hydrogen,3]
Degeneracy_Pressure[Hydrogen,5] = S * Degeneracy_Pressure[Hydrogen,4]
Degeneracy_Pressure[Methane,1] = 0.1
Degeneracy_Pressure[Diamond,2] = SDP/CarbonRadius**P
Degeneracy_Pressure[Diamond,3] = MF * (HDP*3**P) # Ground shell
Degeneracy_Pressure[Diamond,4] = S*S* Degeneracy_Pressure[Diamond,3]
Degeneracy_Pressure[Diamond,5] = S * Degeneracy_Pressure[Diamond,4]
Degeneracy_Pressure[Water,1] = SDP/OxygenRadius**P
Degeneracy_Pressure[Water,2] = MF * (SDP*4**P)/S
Degeneracy_Pressure[Water,3] = MF * (HDP*4**P) # Ground shell
Degeneracy_Pressure[Water,4] = S*S* Degeneracy_Pressure[Water,3]
Degeneracy_Pressure[Water,5] = S * Degeneracy_Pressure[Water,4]
Degeneracy_Pressure[Rock,1] = SDP/OxygenRadius**P
Degeneracy_Pressure[Silicon,2] = MF * (SDP*6.5**P)/S
Degeneracy_Pressure[Silicon,3] = MF * (HDP*6.5**P) # Ground shell
Degeneracy_Pressure[Silicon,4] = S*S* Degeneracy_Pressure[Silicon,3]
Degeneracy_Pressure[Silicon,5] = S * Degeneracy_Pressure[Silicon,4]
Degeneracy_Pressure[Iron,1] = MF * SDP/(IronRadius**P)
Degeneracy_Pressure[Iron,2] = S * Degeneracy_Pressure[Iron,1]
Degeneracy_Pressure[Iron,3] = S * Degeneracy_Pressure[Iron,2]
Degeneracy_Pressure[Iron,4] = MF * (HDP*13**P) # Ground shell
Degeneracy_Pressure[Iron,5] = S*S* Degeneracy_Pressure[Iron,4]
for x in range(1, Elements+1, 1):
for y in range(6, Nl+1, 1):
Degeneracy_Pressure[x,y]=S*Degeneracy_Pressure[x,y-1]
for x in range(0, Elements, 1):
for y in range(1, Nl+1, 1):
Element[x*Nl+y]=x+1
Level[x*Nl+y]=y
# Sorting:
x=1
while x < Elements*Nl+1:
if ElementDensity[Element[x],Level[x]]<ElementDensity[Element[x-1],Level[x-1]]:
Temporary=Element[x]
Element[x] = Element[x-1]
Element[x-1] = Temporary
Temporary=Level[x]
Level[x] = Level[x-1]
Level[x-1] = Temporary
x-=1
else:
x+=1
for L in range(0, Elements*Nl+1, 1):
Layer[Element[L],Level[L]] = L
Density[L]=ElementDensity[Element[L],Level[L]] # change density to earth masses per earth volume
DP[L]=Degeneracy_Pressure[Element[L],Level[L]]
for L in range(0, Elements*Nl+1, 1):
PostCompressionLayer[L]=Layer[Element[L],Level[L]+1]
OriginalLayer[L]=Layer[Element[L],1]
for L in range(1, Elements*Nl+1, 1):
for LL in range(1, Elements*Nl+1, 1):
Ratio[L,LL]=Density[L]/Density[LL]
L=1
print ('\nLayer', 'ElementDensity', 'Element', 'Level', 'Degeneracy-Pressure', 'PostCompressionLayer')
while L < Elements*Nl+1:
print (L, '\t', Density[L], '\t', Element[L], '\t', Level[L], '\t', DP[L], '\t', PostCompressionLayer[L])
L+=1
print ('\n')
fig=plt.gcf()
fig.set_size_inches(14.4,
9.6)
plt.title('Exoplanet structure',
fontsize=16)
plt.ylabel('Jupiter Radii',
color='#000000',
fontsize=16)
plt.xlabel('Jupiter Masses',
color='#000000',
fontsize=16)
#plt.axis([2, 10, 0.5, 2])
#plt.axis([0, 0.02, 0, 0.2])
#plt.axis([0, 0.2, 0, 0.6])
#plt.axis([0.06, 1.5, 0.01, 2]) # log log plot
ax=plt.gca()
ax.set_xscale('log')
ax.set_yscale('log')
plt.minorticks_on()
plt.grid(which='major',
axis='both',
color='black',
linestyle='-')
plt.grid(which='minor',
axis='both',
color='#AAAAAA',
linestyle='dotted')
#plt.gcf().set_facecolor('green')
#ax.set_facecolor('k')
Colors=['#000000',
'#FFFF00', '#FF00FF', '#0000FF', '#00FF00', '#FF8000', '#FF0000',
'#000000']
Colors2=['#000000',
'#FFFFFF', '#000000', '#FFFFFF', '#000000', '#FFFFFF', '#000000', '#FFFFFF', '#000000',
'#FFFFFF', '#000000', '#FFFFFF', '#000000', '#FFFFFF', '#000000', '#FFFFFF', '#000000',
'#FFFFFF', '#000000', '#FFFFFF', '#000000', '#FFFFFF', '#000000', '#FFFFFF', '#000000',
'#000000']
Marker=['o',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o', 'o', 'o', 'o',
'+', '+', '+', '+',
'o']
# **********************************************************************
# **********************************************************************
# **********************************************************************
# **********************************************************************
# **********************************************************************
# **********************************************************************
# **********************************************************************
# **********************************************************************
def GBE(InnerMass, OuterMass, GhostMass, InnerRadius, OuterRadius):
if InnerRadius>0:
InnerMass -= GhostMass
OuterMass += GhostMass
return ((InnerMass/InnerRadius - InnerMass/OuterRadius) + (0.5*OuterMass/OuterRadius - (0.5*OuterMass*InnerRadius*InnerRadius)/(OuterRadius**3)))
else:
return 0.5*OuterMass/OuterRadius
StartTime=time.time()
Time=0
Multiplier=1.05**1
Mal=1
CoreLevel=1
Collapse1=0
Collapse2=0
T=1
if InitialRockVolume==0 and InitialIronVolume==0:
for L in range(0, Elements*Nl+1, 1):
Volume[L]=0.1*Delta2[1,Element[L],Level[L]] # our initial volumes
Delta[L]=Volume[L] # because we added this much to zero
if Volume[L]>0:
Finished[L]=0
else:
Finished[L]=1
else:
Volume[Layer[Rock,1]]=InitialRockVolume
Delta[Layer[Rock,1]]=InitialRockVolume
Volume[Layer[Iron,1]]=InitialIronVolume
Delta[Layer[Iron,1]]=InitialIronVolume
MaxTime=330
Time=1
while Time < (MaxTime+1):
StillCollapsing=1
while StillCollapsing:
StillCollapsing=0
Compressing=1
while Compressing:
Compressing=0
Count[0]+=1
if Count[0]>1000000:
print ('Taking too long')
print ('level\tfinished\telement\tlevel\tcount\tVolume')
for L in range(1,Elements*Nl+1,1):
print (L, Finished[L], Element[L], Level[L], Count[L], Volume[L])
plt.savefig('Exoplanet.Structure.All.png', dpi=100) # must come before the show command. Width in pixels = dpi * 6.4
duration = 1000 # milliseconds
freq = 1000 # hz
winsound.Beep(freq,duration)
raise SystemExit()
for Lyr in range(Elements*Nl, 0, -1): # bottom to top
TotalVolume[Lyr]=TotalVolume[Lyr+1] + Volume[Lyr]
Radius[Lyr]= TotalVolume[Lyr]**0.333333
GhostMass[Lyr]= TotalVolume[Lyr+1]*Density[Lyr]
Mass[Lyr]= Density[Lyr]*Volume[Lyr]
TotalMass[Lyr]= Mass[Lyr]+TotalMass[Lyr+1]
if Volume[Lyr]>0:
# 1 g's * 1 earth radius * 1 earth mass per earth volume in bar = 3,450,000 bar
Pressure[Lyr]=Density[Lyr]*GBE(TotalMass[Lyr+1], Mass[Lyr], GhostMass[Lyr], Radius[Lyr+1], Radius[Lyr])
else:
Pressure[Lyr]=0
# BindingEnergy[0]=(TotalMass[0]/Radius[1])
# TotalBindingEnergy[0]=BindingEnergy[0]
# for x in range(1, Nl*Elements+1, 1):
# BindingEnergy[x]=GBE(TotalMass[x+1], Mass[x], GhostMass[x+1], Radius[x+1], Radius[x])
# TotalBindingEnergy[x]=BindingEnergy[x]+TotalBindingEnergy[x-1]
for L in range(1, Elements*(Nl-1)+1, 1): # top to bottom
TF=1
if Element[L]==5 and Level[L]==1 and CoreTemperature[CoreLevel]<CuriePoint[Iron,1]:
TF=MMF
if Volume[L]>0 and (Pressure[L]+TotalPressure[L-1])>DP[L]*TF:
TotalPressure[L]=DP[L]*TF+TotalPressure[L-1]
else:
TotalPressure[L]=Pressure[L]+TotalPressure[L-1]
# if Element[L]==3 and Level[L]==3 and CoreTemperature[CoreLevel]<CuriePoint[Water,3]:
# TF=MMF
# if Element[L]==1 and Level[L]>1:
# T=MMF
if Volume[L]>0 and (Pressure[L]+TotalPressure[L-1])>DP[L]*TF:
PCL=PostCompressionLayer[L]
OL=OriginalLayer[L]
Count[L]+=1
Decrease=0.001*Delta[OL]*Ratio[OL,L]
if Decrease>Volume[L]:
Decrease=Volume[L]
if Element[L]==Rock and Level[L]==1:
Volume[L]=Volume[L]-Decrease
Volume[Layer[Water,2]]+=Decrease*0.5*Ratio[L,Layer[Water,2]]
Volume[PCL]+=Decrease*0.5*Ratio[L,PCL]
else:
Volume[L]-=Decrease
Volume[PCL]+=Decrease*Ratio[L,PCL]
Compressing=1
if Collapse1==0 and Volume[Layer[Iron,2]]>0:
Percent=1
Volume[Layer[Iron,3]]=Volume[Layer[Iron,3]]+Percent*Volume[Layer[Iron,2]]*Ratio[Layer[Iron,2],Layer[Iron,3]]
Volume[Layer[Iron,2]]=(1-Percent)*Volume[Layer[Iron,2]]
Collapse1=1
StillCollapsing=1
InitialCompressedIronVolume=Volume[Layer[Iron,3]]
print ('collapse1')
if Collapse1==1 and Collapse2==0 and Volume[Layer[Water,4]]>0:#Volume[Layer[Water,2]]>0:
print ('collapse2')
# Volume[Layer[Iron,4]]=Volume[Layer[Iron,4]]+Volume[Layer[Iron,3]]*Ratio[Layer[Iron,3],Layer[Iron,4]]
# Volume[Layer[Iron,3]]=0
# Volume[Layer[Iron,4]]=Volume[Layer[Iron,4]]+Volume[Layer[Iron,2]]*Ratio[Layer[Iron,2],Layer[Iron,4]]
# Volume[Layer[Iron,2]]=0
# Volume[Layer[Iron,4]]=Volume[Layer[Iron,4]]+Volume[Layer[Iron,1]]*Ratio[Layer[Iron,1],Layer[Iron,4]]
# Volume[Layer[Iron,1]]=0
# Volume[Layer[Rock,4]]=Volume[Layer[Rock,4]]+Volume[Layer[Rock,3]]*Ratio[Layer[Rock,3],Layer[Rock,4]]
# Volume[Layer[Rock,3]]=0
# Volume[Layer[Rock,4]]=Volume[Layer[Rock,4]]+Volume[Layer[Rock,2]]*Ratio[Layer[Rock,2],Layer[Rock,4]]
# Volume[Layer[Rock,2]]=0
# Volume[Layer[Rock,4]]=Volume[Layer[Rock,4]]+Volume[Layer[Rock,1]]*Ratio[Layer[Rock,1],Layer[Rock,4]]
# Volume[Layer[Rock,1]]=0
# Volume[Layer[Water,4]]=Volume[Layer[Water,4]]+Volume[Layer[Water,3]]*Ratio[Layer[Water,3],Layer[Water,4]]
# Volume[Layer[Water,3]]=0
# Volume[Layer[Water,4]]=Volume[Layer[Water,4]]+Volume[Layer[Water,2]]*Ratio[Layer[Water,2],Layer[Water,4]]
# Volume[Layer[Water,2]]=0
# Volume[Layer[Water,4]]=Volume[Layer[Water,4]]+Volume[Layer[Water,1]]*Ratio[Layer[Water,1],Layer[Water,4]]
# Volume[Layer[Water,1]]=0
# Volume[Layer[Methane,4]]=Volume[Layer[Methane,4]]+Volume[Layer[Methane,3]]*Ratio[Layer[Methane,3],Layer[Methane,4]]
# Volume[Layer[Methane,3]]=0
# Volume[Layer[Methane,4]]=Volume[Layer[Methane,4]]+Volume[Layer[Methane,2]]*Ratio[Layer[Methane,2],Layer[Methane,4]]
# Volume[Layer[Methane,2]]=0
# Volume[Layer[Methane,4]]=Volume[Layer[Methane,4]]+Volume[Layer[Methane,1]]*Ratio[Layer[Methane,1],Layer[Methane,4]]
# Volume[Layer[Methane,1]]=0
#
## Volume[Layer[Methane,5]]=Volume[Layer[Methane,5]]+Volume[Layer[Methane,4]]*Ratio[Layer[Methane,4],Layer[Methane,5]]
## Volume[Layer[Methane,4]]=0
#
## Volume[Layer[Hydrogen,3]]=Volume[Layer[Hydrogen,3]]+Volume[Layer[Hydrogen,1]]*Ratio[Layer[Hydrogen,1],Layer[Hydrogen,3]]
## Volume[Layer[Hydrogen,1]]=0
Percent=1
Volume[Layer[Water,5]]=Volume[Layer[Water,5]]+Percent*Volume[Layer[Water,4]]*Ratio[Layer[Water,4],Layer[Water,5]]
Volume[Layer[Water,4]]=(1-Percent)*Volume[Layer[Water,4]]
Percent=1
Volume[Layer[Methane,5]]=Volume[Layer[Methane,5]]+Percent*Volume[Layer[Methane,4]]*Ratio[Layer[Methane,4],Layer[Methane,5]]
Volume[Layer[Methane,4]]=(1-Percent)*Volume[Layer[Methane,4]]
Collapse2=1
# StillCollapsing=1
Mal=5
for L in range(1, Nl*Elements+1, 1): # top to bottom
if Volume[L]>0:
plt.plot(TotalMass[1]/JM,
Radius[L]/JR,
marker='o',
color=Colors[Element[L]],
markerfacecolor=Colors2[Level[L]],
markersize=3,
markeredgewidth=1,
markeredgecolor=Colors[Element[L]],
alpha=1.0,
fillstyle='full')
print (Time, '\t', end="")
Time+=1
if Time<(MaxTime+1):
# if Mal==4 and TotalMass[1]>Delta2[5, 0, 0]:
# Mal=5
if Mal==3 and TotalMass[1]>Delta2[4, 0, 0]:
Mal=4
elif Mal==2 and TotalMass[1]>Delta2[3, 0, 0]:
Mal=3
elif Mal==1 and TotalMass[1]>Delta2[2, 0, 0]:
Mal=2
AddedMass=0
TargetMass=TotalMass[1]*Multiplier
for x in range(1, Elements+1, 1):
for y in range(1, Nl+1, 1):
AddedMass=AddedMass+Delta2[Mal,x,y]*ElementDensity[x,y]
Increase=(TargetMass-TotalMass[1])/AddedMass
for L in range(1, Elements*Nl+1, 1): # top to bottom
Delta[L]=Delta2[Mal,Element[L],Level[L]]*Increase
Volume[L] = Volume[L] + Delta[L]
if Delta[L]>0:
Finished[L]=0
x=Nl
while x>0:
if Volume[Layer[Iron,x]]>0:
CoreLevel=x
x=0
x-=1
print (Delta[1])
plt.annotate('Iron below\nCurie point',
ha='center',
xy=(0.0008, 0.04),
xytext=(0.0003, 0.015),
arrowprops=dict(color='#000000',
shrink=0.1,
width=1,
headwidth=6,
headlength=6))
plt.annotate('Iron above Curie point',
ha='center',
xy=(0.005, 0.02),
xytext=(0.05, 0.02),
arrowprops=dict(color='#000000',
shrink=0.1,
width=1,
headwidth=6,
headlength=6))
plt.annotate('Hydrogen fusion\nbegins',
ha='center',
xy=(80, 1.5),
xytext=(80, 2),
arrowprops=dict(color='#000000',
shrink=0.1,
width=1,
headwidth=6,
headlength=6))
plt.annotate('helium fusion\nbegins',
ha='center',
xy=(700, 1.5),
xytext=(700, 2),
arrowprops=dict(color='#000000',
shrink=0.1,
width=1,
headwidth=6,
headlength=6))
plt.annotate('Diamond and water floating above metallic hydrogen',
xy=(5, 0.9),
xytext=(10, 0.4),
arrowprops=dict(color='#000000',
shrink=0.1,
width=1,
headwidth=6,
headlength=6))
plt.annotate('Pure helium atmosphere',
ha='center',
xy=(3, 1.5),
xytext=(3, 2.1),
arrowprops=dict(color='#000000',
shrink=0.1,
width=1,
headwidth=6,
headlength=6))
plt.text(0.002,
0.4,
'Electron shell\ndegeneracy pressure\n' + r'$P = \frac{n}{r^6}$' + '\nn = number of electrons',
horizontalalignment='center',
fontsize=16 )
plt.plot(1/JM,
0.017,
marker='o',
color='#0000FF',
markerfacecolor='#0000FF',
markersize=3,
markeredgewidth=1,
markeredgecolor='#0000FF',
alpha=1.0,
fillstyle='full')
plt.plot(1/JM,
0.04876,
marker='o',
color='#0000FF',
markerfacecolor='#0000FF',
markersize=3,
markeredgewidth=1,
markeredgecolor='#0000FF',
alpha=1.0,
fillstyle='full')
#
#print ('level\tfinish\telement\tlevel\tcount\tdelta')
#for L in range(1,Elements*Nl+1,1):
# print (L, '\t', Finished[L], '\t', Element[L], '\t', Level[L], '\t', Count[L], Delta[L])
EndTime=time.time()
print ('\nNumber of loops:', Count[0])
print ('Time:', EndTime-StartTime)
print ('\nLayer', 'Count', 'Element', 'Level', 'Volume', 'TotalPressure')
L=1
while L < Elements*Nl+1:
if Volume[L]>0:
print (L, '\t', Count[L], '\t', Element[L], '\t', Level[L], '\t', Volume[L], '\t', TotalPressure[L])
L+=1
plt.savefig('Exoplanet.Structure.All.png',
dpi=100) # must come before the show command. Width in pixels = dpi * 6.4
duration = 1000 # milliseconds
freq = 440 # hz
winsound.Beep(freq,duration)
plt.show()
