diff --git a/Time-Stepping Naked Planet Model/Grading/GWmod1.py b/Time-Stepping Naked Planet Model/Grading/GWmod1.py new file mode 100644 index 0000000..90356e0 --- /dev/null +++ b/Time-Stepping Naked Planet Model/Grading/GWmod1.py @@ -0,0 +1,76 @@ +#My first tijme making a plot for climate modeling. +#Hopefully, this works. +#Note that using spreadsheets for this course is optional; +#the only way to get credit is with Python3! + +"""The background for this material can be found in Sections 2 and 3 of Part I +of this class. Joules are units of energy; Watts are energy flow (J/s). +The temperature of a planet is determined by balancing energy fluxes in and out of a planet. +Incoming solar heat is determined by L*(1-albedo)/4, +and outgoing infrared is calculated as epsilon * sigma * T^4.""" + +"""The heat capacity of the surface planet depends on the water thickness. +Using the fact that 1 gram of water warms by 1 K with a 4.2J (1 cal) of heat, +and the density of water, to calculate how many Joules of energy it takes +to raise the temperature of 1 m2 of the surface, given some water depth +in meters.""" + +"""The goal is to numerically simulate how the planetary temperature of a +naked planet would change through time as it approaches equilibrium +(the state at which it stops changing, which we calculated before). +The planet starts with some initial temperature. +The “heat capacity” (units of Joules / m2 K) of the planet is set by +a layer of water which absorbs heat and changes its temperature. +If the layer is very thick, it takes a lot more heat (Joules) to change the temperature. +The differential equation you are going to solve is:""" +#d(HeatContent)/dt = L*(1-alpha)/4 - epsilon * sigma * T^4 +"where the heat content is related to the temperature by the heat capacity:" +#T[K] = HeatContent [J/m2] / HeatCapacity [J/m2 K] +"""The numerical method is to take time steps, extrapolating the heat content +from one step to the next using the incoming and outgoing heat fluxes, +same as you would balance a bank account by adding all the income and subtracting all the expenditures +over some time interval like a month. +The heat content "jumps" from the value at the beginning of the time step, +to the value at the end, by following the equation:""" +#HeatContent(t+1) = HeatContent(t) + dHeatContent/dT * TimeStep +"""This scheme only works if the time step is short enough so that +nothing too huge happens through the course of the time step. +Set the model up with the following constants:""" +timeStep = 40 # years; convert to seconds. +time = 0 # years; this shall denote how much time has passed in total. +T = 0 # Kelvin, according to question demands. +waterDepth = 2000 # meters; multiply by cH2O and water density. +L = 1350 # W m-2, same as J s-1 m-2 (Solar Constant) -> multiply by no. of seconds in each timestep +albedo = 0.3 # alpha +epsilon = 1 +sigma = 5.67E-8 # W m-2 K-4 +cH2O = 4200 # specific heat capacity of water (units of 4200 J kg-1 K-1). +rhoH2O = 1000 # Density is 1000 kg m-3. +energy_per_m2 = waterDepth * cH2O * rhoH2O #in J m-2 K-1 +#heat influx and outflux must be calculated along the way. +import numpy +import matplotlib.pyplot +#Getting our lists needed with the very first variables already in: +nSteps = input("Please use 20 or more") +#apparently the autograder that the MOOC uses also mistakes strings for operators. +nSteps = int(nSteps) +time_list = [time] +temperature_list = [T] +HeatContent = 0 #only true when T = 0 Kelvin, this is the heat content of the Earth. +difference = 0 +outflux = 0 +for x in range(nSteps): + time = time + timeStep + time_list.append(time) + difference = (L * (1 - albedo) / 4 - (epsilon * sigma * T**4)) #W m-2 + HeatContent = HeatContent + difference * 365 * 24 * 60 * 60 * timeStep #J m-2 + T = (HeatContent / energy_per_m2) #Kelvin + outflux = epsilon * sigma * T**4 # as per instructions + temperature_list.append(T) +#derive answers and end with: +matplotlib.pyplot.plot(time_list, temperature_list) +matplotlib.pyplot.show() +#Give the last temperature derived by the model: +#print(temperature_list[-1], outflux) +#The tolerances in the grader for temperatures are 2-3 degrees C, +#and for the heat fluxes, +- 0.1 W/m2. diff --git a/Time-Stepping Naked Planet Model/Grading/Naked_Earth.py b/Time-Stepping Naked Planet Model/Grading/Naked_Earth.py new file mode 100644 index 0000000..8a63aea --- /dev/null +++ b/Time-Stepping Naked Planet Model/Grading/Naked_Earth.py @@ -0,0 +1,32 @@ +import matplotlib.pyplot as plt + +timeStep = 10 # years +waterDepth = 4000 # meters +L = 1350 # Watts/m2 +albedo = 0.3 +epsilon = 1 +sigma = 5.67E-8 # W/m2 K4 +#nSteps = int(input("")) +nSteps=200 + +heat_cap=waterDepth*4.2E6 #J/Kg/deg +time_list=[0] +temperature_list=[400] +heat_content=heat_cap*temperature_list[0] +heat_in=L*(1-albedo)/4 +heat_out=0 +secs_in_yr=3.14E7 + +for t in range(0,nSteps): + time_list.append(timeStep+time_list[-1]) + heat_out=epsilon*sigma*pow(temperature_list[-1],4) + heat_content=heat_content+(heat_in-heat_out)*timeStep*secs_in_yr + temperature_list.append(heat_content/heat_cap) + heat_out = epsilon * sigma * pow(temperature_list[-1], 4) +print(temperature_list[-1]) +#print(temperature_list[-1], heat_out) +plt.plot(time_list, temperature_list, marker='.', color="k") +plt.xlabel("Time (years)") +plt.ylabel("Temperature (degrees Kelvin)") +plt.show() + diff --git a/Time-Stepping Naked Planet Model/Grading/naked_planet_tosubmit.py b/Time-Stepping Naked Planet Model/Grading/naked_planet_tosubmit.py new file mode 100644 index 0000000..4d567f2 --- /dev/null +++ b/Time-Stepping Naked Planet Model/Grading/naked_planet_tosubmit.py @@ -0,0 +1,68 @@ +# %% + + +# %% +import numpy as np +import matplotlib.pyplot as plt + +# Set up constants as given: +timeStep_years = 100 # years +waterDepth = 4000 # meters +L = 1350 # W/m² +albedo = 0.3 +epsilon = 1 +sigma = 5.67e-8 # W/m² K^4 + +# The incoming solar flux is constant: +incoming_flux = L * (1 - albedo) / 4 # W/m² + +# Calculate the heat capacity (J/m² K) +# (water density = 1000 kg/m³, specific heat = 4184 J/(kg·K)) +heat_capacity = waterDepth * 1000 * 4184 + +# Read the number of time steps from stdin +nSteps = int(input("")) + +# Set the initial conditions. +# (Initial temperature is chosen in Kelvins; 0 K is allowed though not physical.) +initial_temperature = 0.0 # Kelvin +heatContent = initial_temperature * heat_capacity # in J/m² + +# Prepare lists for time and temperature (for plotting, if desired) +time_list = [0] +temperature_list = [initial_temperature] + +# Convert the time step to seconds. +dt = timeStep_years * 365 * 24 * 3600 + +# Integration loop: update heat content and temperature each step. +for i in range(nSteps): + # Current temperature in Kelvin: + T = heatContent / heat_capacity + + # Differential equation: dHeat/dt = incoming_flux - ε·σ·(T)⁴ + dH_dt = incoming_flux - epsilon * sigma * (T**4) + + # Update the heat content over the time step: + heatContent = heatContent + dH_dt * dt + + # Append new time and temperature to lists: + time_list.append((i+1) * timeStep_years) + temperature_list.append(heatContent / heat_capacity) + +# Compute the final temperature (in Kelvin) and the corresponding outgoing heat flux. +final_temperature = heatContent / heat_capacity +final_heat_flux = epsilon * sigma * (final_temperature**4) + +# Output only the final temperature and heat flux. +print(final_temperature, final_heat_flux) + +# Plotting code +plt.plot(time_list, temperature_list) +plt.show() + + +# %% + + + diff --git a/Time-Stepping Naked Planet Model/naked_planet.py b/Time-Stepping Naked Planet Model/naked_planet.py new file mode 100644 index 0000000..66cce23 --- /dev/null +++ b/Time-Stepping Naked Planet Model/naked_planet.py @@ -0,0 +1,51 @@ +import numpy +import matplotlib.pyplot + + +# Inputs +timeDuration = 2000 # years +timeStep = 10 # years +waterDepth = 4000 # meters +initialTemp = 400 + +# Constants +secondToYear = 60*60*24*365 +specificWaterHeatCapacity = 4.2E3 # kg/m^3 +waterDensity = 1000 # J/kg*m^3 +L = 1350 # Watts/m2 +albedo = 0.3 +epsilon = 1 +sigma = 5.67E-8 # W/m2 K4 + + +# Calc amount of steps for given duration and resolution +nSteps = int(timeDuration/timeStep) + +initialHeat = initialTemp * (waterDepth*waterDensity*specificWaterHeatCapacity) +initialRadiatedHeat = epsilon * sigma * pow(initialTemp,4) +T = [initialTemp] # K +Q = [initialHeat] # Watt/m^2 +time = [0] # years +heatOut = [initialRadiatedHeat] # Watt/m^2 +# print(initialRadiatedHeat) + +for i in range(nSteps): + # Calculate the heat change rate + dQ_dt = (L * (1 - albedo)/4) - epsilon * sigma * pow(T[-1],4) + # Add heat change over step period to total heat in system + Q.append(Q[-1] + dQ_dt * timeStep * secondToYear) + # Calculate temperature of planet based on heat in system + T.append(Q[-1] / (waterDepth*waterDensity*specificWaterHeatCapacity)) + # Calculate only the outgoing heat + heatOut.append(epsilon * sigma * pow(T[-1],4)) + # Save timestep for plotting + time.append(time[-1] + timeStep) + +print(T[-1]) + + +matplotlib.pyplot.plot(time, T) +matplotlib.pyplot.title("Time-dependent temperature of a naked planet") +matplotlib.pyplot.xlabel('Time (years)') +matplotlib.pyplot.ylabel('Temperature (K)') +matplotlib.pyplot.show() \ No newline at end of file diff --git a/fortran setup/main.f95 b/fortran setup/main.f95 deleted file mode 100644 index 1be38f6..0000000 --- a/fortran setup/main.f95 +++ /dev/null @@ -1,4 +0,0 @@ -program main - implicit none - write(*,*) "Hello world!" -end program main