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