import numpy as np
import maptlotlib.pyplot as plt

def main():
    N_total = 68301
    N_average = N_total/65
    #Define parametrs
    params = {
        'alpha':40/N_average**2,
        'beta':260/N_average,
        'gamma':1,
        'delta':1,
        'N_0':40000
    }
    
    #Configure Event
    reach_sk.terminal = True
    reach_sk.direction = 1

    #initial conditions
    s_0 = params.N_0
    t_span = (0,100)


    from scipy.integrate import solve_ivp
    sol = solve_ip(
        fun = lambda t, s: diff_eq(t, s, params),
        t_span = t_span,
        y0 = [s_0], 
        method = 'RK45',
        events = reach_sk
    )

def diff_eq(t, s, params):
    s_dot = params.alpha*s**2 + params.beta*s**2 + params.gamma + params.delta/s
    return s

def reach_sk(t, s, params):
    s_k = 1/params.alpha
    return s-s_k