##############################
#                            #
#  PP04 ice age simulator    #
#                            #
##############################

######################################
#
#  V:  Overall sea ice
#  A:  Antarctic land ice
#  C:  Atmospheric CO2
#  t:  Time in kilo years
#
######################################

import numpy as np
import scipy as sci
import matplotlib.pyplot as plt
from scipy.integrate import odeint

#####################################
#
#  Forcing insolation parameters
#
#####################################

########################################################
#
#  Insolation forcing is  mu cos(om*t)
#  
#  Representative physical values are:
#
#  om = 0.13
#  mu = 0.4
#
#########################################################

omm = input(" Give the value of omega " )
muu = input(" Give the value of mu ")

om = float(omm)
mu = float(muu)


############################################################################
############################################################################


def ff(x,t,om,mu):
##################

########################################
#                                      #
#   RHS of the PP04 model              #
#                                      #
########################################



    L   = np.array([[-1/15, 0, -1.3/15],[1/12, -1/12, 0],[-0.5/5, 0, -1/5]])
    e   = np.array([-0.5/15, 0, 0.15/5])
    bp  = np.array([0.8/15, 0, 0.4/5])
    bm  = np.array([0.8/15, 0, 1.1/5])

    


###########################################################################
#
#  Test for the region of glacial: (F > 0)  vs. inter-glacial: (F < 0)
#
###########################################################################

    dd  = 0.27
    cc  = np.array([0.3, -0.7, 0])
    F   = np.dot(cc,x) + dd 

###################################################
#  F measures the amount of CO2 in the oceans
###################################################
#
# ‘Discontinuity’ function
#
########################################

    eta = 2000     
    H = 0.5*(1 + np.tanh(eta*F))

#####################################################
#
#  ODE
#
#####################################################

    f = L.dot(x) + bm*(1-H)+bp*H + e*mu*np.sin(om*t)

    
    return f

######################################################
######################################################


######################################################
#
#  Intial values (these are not very important)
#
######################################################

u0 = np.array([0.5,0.5,0.5])

tv = np.linspace(0,400,400)


######################################################
#
#  Integrate the PP04 model
#
######################################################


u  = odeint(ff,u0,tv,args=(om,mu))


######################################################


######################################################
#
#   Plot the solution
#
######################################################


plt.plot(tv,u)

plt.xlabel('time')
plt.ylabel('[V,A,C]')
plt.title('PP04 Model of the ice ages')

plt.show()