0

First let me apologize if this is a simple question, or it has been answered before. The truth is that I don't even know how to run a search about my question.

Let's say I have the following set of coupled ODE:

dN/dt = (hN(t) + iQ(t))*G(t)

dQ/dt = (jN(t) + lQ(t))*G(t)

This is the code i used to plot the ODEs:

import numpy 
import matplotlib.pyplot as plt
from scipy import constants
from scipy.integrate import odeint
plt.ion()

#-------------------------------------INITIAL CONSTANT PARAMETERS-------------------------------------------

wavelenght      = 1908e-9 #meter
wp0             = 300e-6 #laser beam diameter
conc            = 1.0 #dopant concentration
crosssection    = 6.5e-14 #spectroscopic absorption cross section in m2
lifetime        = 230e-6 #upper lifetime
cavity_l        = 200e-3 #total lenght of cavity in mm
gain_l          = 10e-3 #lenght of gain medium
n_gain          = 1.82 #refraction index gain medium
R_OC            = 0.3 #reflectivity OC
additional_loss = 0.05 #additional losses
loss_time_max   = 10.0e5 #max loss instroduced by q-switch
r               = 10 #times above threshold laser is pumped before q-switch opens
t0              = 0
tf              = 20e-12
tpulse          = 1e-9
ttotal          = 10e-9

#-------------------------------------INITIAL VARIABLE PARAMETERS-------------------------------------------

#number density from concentration percentage
Ntotalyag             = conc*3*4.55/((3*88.9 + 5*27.0 + 12*16.0)*constants.m_p)
#nominator: 3 at of Y * mass density of Y3Al5O12
#denominator: mass of the Y3Al5O12 unit, calculated from the relative atomic weights and the proton mass

beam_area             = numpy.pi*(wp0/2)*(wp0/2)

roundtrip_time        = 2.0*((cavity_l-gain_l)/constants.c)+(n_gain*gain_l/constants.c) #time for light to travel back and forth inside cavity

losscoef              = - numpy.log(R_OC) + additional_loss 

popinversionthreshold = losscoef * beam_area/(2.0 * crosssection)

Wp = 0.0 #pumping rate

#-------------------------------------FUNCTIONS-------------------------------------------

def f(y, t, L):
    q      = y[0]
    DeltaN = y[1]
    if t > tf:
        losscoef_t = 0.0
    else:
        losscoef_t = loss_time_max - ((loss_time_max - losscoef)/(tf))*t
    # gupta Handbook of photonics ch about micro laser
    L= L.append(losscoef_t)
    f0 = q*((DeltaN/popinversionthreshold)-1.0) * losscoef_t/roundtrip_time
    f1 = Wp*(Ntotalyag-DeltaN) - q*((DeltaN*losscoef_t)/(popinversionthreshold*roundtrip_time)) - DeltaN/lifetime
    return [f0, f1]

#-------------------------------------ODE SOLVING 1-------------------------------------------
# initial conditions
q0 = 0.01               # initial photon population
DeltaN0 = r * popinversionthreshold   # initial inverted population
y0 = [q0, DeltaN0]       # initial condition vector
t  = numpy.linspace(0, ttotal, 1000)   # time grid

L = []

# solve the DEs
soln = odeint(f, y0, t, args=(L,))
Q    = soln[:, 0] 
N    = soln[:, 1] 

#-------------------------------------PLOTS---------------------------------------------------

plt.figure()
plt.plot(t, Q, label='Photon')
plt.plot(t, N, label='Population')
plt.legend(loc=0)

The above does not work, since i assumed i could calculate G(t) for each timestep t, but i noticed len(L) is not the same as Q or N. Basically my question is, how to i calculate a parameter on the ODE that varies with the timestep t?

Thanks a lot

Hamilton Pupo
  • 1
  • 1. Not sure why you use L = L.append(x) rather than just L.append(x). L is not the same length of Q or N because the function will be called many times between the timesteps you requested solutions at. One could, presumably, have a second list T, where each time you append a value to L you also append the time. You will then have L vs T for all the times the function was evaluated. (Yes, ugly). – Jon Custer Jun 19 '15 at 19:05
  • thanks for the feedback. Actually right after i post this i did a bit more research and realized my whole code is messed up. I will be closing this thread as shortly, because the entire post does not make too much sense. – Hamilton Pupo Jun 22 '15 at 13:38
  • Somehow, for me, the act of reducing code to something clean and simple to post makes whatever problem I'm having become clear and obvious. Glad to have been of some help (perhaps). Good luck! – Jon Custer Jun 22 '15 at 13:42

0 Answers0