####  本算法来源于lichagnsheng老师最新的书，结果看上去似乎偏小，需要进一步的确定
%matplotlib inline
import numpy as np

import matplotlib.pyplot as plt
import unittest
import random
TCfilm=0.25      #film thermal conductivity, 0.25 J/(s*m*K)
Dfilm=0.001      #film thickness, 0.001 m //0.0075 mm
dt = 86400      #time, second	
# 1 calorie = 4.184 J
SHsoil = 0.2    # soil heat capacity, 0.2 calorie/g/degree C

adjT = np.random.random(100)/2*100
#Surf_Temp = np.random.random(100)/2*100

#Surf_Temp = adjT-20.0
Surf_Temp = np.random.random(100)/2*100

# assume 2 cm per layer
# soil mass = soil depth * soil density * area in ha... 1.3(g/cm3) * 2(cm) * 100000000(cm2) / 1000
soil_depth = 0.02 # meter
soil_bulk_density = 1.3 #(g/cm3) /1300 kg/m3

Soil_mmm = 0.02 * 10000 * 1300  # kg

print Soil_mmm
#Soil_mmm = 114660   # soil mass kg/ha// 

Msoil = 220200        #

SHsoil = 0.0478 * Soil_mmm
dH = (adjT - Surf_Temp) * Dfilm * TCfilm * dt

dTTT = dH * SHsoil


#dTTT = (dHHH / 4.184*1000) / (SHsoil * Soil_mmm)

#plt.plot(adjT,label="air temperature")
#plt.plot(Surf_Temp,label="Previous day surface temperature/air tmpe+20")
plt.plot(dTTT,label="Changed temperature")
plt.legend(loc='best')
plt.show()

260000.0