File talk:Planetary interior structure.png
Latest comment: 5 years ago by Just granpa
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()
