# -*- coding: utf-8 -*-
"""
Created on Mon Feb 17 11:07:09 2020

@author: neil
"""

# Numerical Modeling in the Geosciences,  Humphrey 2022
# 2nd problem, include radioactive heat generation, 
# ***************************************************************************
# Problem: steady state temperature field in a 20000 m thick layer of rock,
# with a surface temperature of 0degrees and a heat flux of 65 mW per square meter
# This includes reasonable radioactive heating
# The upper BC eqn is q = -K dT/dz
# ***************************************************************************

import numpy as np
import matplotlib.pyplot as plt

K = 3               # thermal conductivity 'granite',
n = 9   			# number of internal nodes, 
thickness = 20000	# thickness of rock layer in vertical, in meters
Rq= .0000015        # nominal radioactive heat generation in upper crust
upperq  = 6.5e-2	# W/m2 upper typical geothermal heat flux (this is a gradient condition)
upperBC = 0 		# temp BCs, in degrees
T = np.zeros([n])
z = np.linspace(0,thickness,len(T))
delx = z[1]-z[0]

# actual program
T[0] = upperBC
T[1] = T[0] + (upperq)*delx/K         # apply the gradient condition by calc the diff between T[0] and T[1]

for j in np.arange(1,n-1):            # integrate from the surface down
    T[j+1] = 2*T[j] - T[j-1] - Rq*delx*delx/K

fig = plt.figure()                    #create a place to plot
ax  = fig.add_subplot(111)            #create axes. 
ax.plot(T, z,'r-o');
ax.invert_yaxis()
ax.set_title('Steady State temperature in radiogenic crust');
ax.set_ylabel('Depth from surface in meters')
ax.set_xlabel('Temperature in degrees')
ax.grid()

plt.show()