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() 

Just granpa (talk) 18:45, 26 October 2018 (UTC)Reply

Return to the file "Planetary interior structure.png".