In [1]:

    

import sympy as sym
import numpy as np
sym.init_printing(use_latex='mathjax')

wlP, wlL, wL, wA = sym.symbols('lambda_P lambda_L w_L w_A', real=True)
TR, tauP, tauL = sym.symbols('T_R tau_P tau_L', real=True)
Toc, etaP = sym.symbols('T_{oc} eta_P', real=True)
DR = sym.Symbol('\Delta R', real=True)
DRns = sym.Symbol('\Delta R_{ns}', real=True)
FsatA, FsatL = sym.symbols('F_{satA} F_{satL}', real=True)
Pstmax, gstmax = sym.symbols('P_{stmax} g_{stmax}', real=True)

# Laser parameters:
v = {wlP: 808e-9,
     Toc: 8.7e-2,
     wlL: 1064e-9,
     wL: 80e-6,
     wA: 140e-6,
     TR: 1 / 100.25e6,
     tauP: 10.6e-12,
     tauL: 100e-6,
     FsatL: 82e-3 / 1e-4,
     FsatA: 60e-6 / 1e-4,
     DR: 1.7e-2,
     DRns: 1.3e-2,
     Pstmax: 4.6}


    

In [2]:

    

import numpy as np
from scipy.optimize import brentq
sym.init_printing(use_latex='mathjax')

def qP(EP, DR, EsatA):
    """Loss at SESAM"""
    S = EP / EsatA
    if type(EP) is float:
        exp = np.exp
    else:
        exp = sym.exp
    return(DR / S * (1 - exp(-S)))

EsatA, EsatL = sym.symbols('E_{satA} E_{satL}', real=True)
Pst, gst, l, EP = sym.symbols('P_{st} g_{st} l E_P', real=True)
PP = sym.symbols('P_P')

# Derived:
v[EsatL] = np.pi * v[wL]**2 * v[FsatL]
v[EsatA] = np.pi * v[wA]**2 * v[FsatA]
v[etaP] = v[wlP] / v[wlL]
v[l] = v[Toc] + v[DRns]
v[gstmax] = v[l] + qP(v[Pstmax] * v[TR], v[DR], v[EsatA])

# Steady-state equations:
gst_ = sym.Eq(gst, l + qP(Pst * TR, DR, EsatA))
Pst_ = sym.Eq(Pst, -EsatL / tauL + etaP * PP / gst)

Pst_.subs(gst_.lhs, gst_.rhs)


    

    
Out[2]:





$$P_{{st}} = - \frac{E_{{satL}}}{\tau_{L}} + \frac{P_{P} \eta_{P}}{\frac{E_{{satA}} \Delta R}{P_{{st}} T_{R}} \left(1 - e^{- \frac{P_{{st}} T_{R}}{E_{{satA}}}}\right) + l}$$

In [3]:

    

# The root of this function is the steady-state intracavity power:
rootfunc = Pst_.subs(gst_.lhs, gst_.rhs).rhs - Pst
rootfunc = rootfunc.subs(v)

# The pump-threshold is:
PPthreshold = (EsatL / tauL * (l + DR) / etaP).subs(v)

N = 100
PPrange = np.linspace(PPthreshold, (Pstmax * Toc / etaP).subs(v), N)
Pstrange = np.zeros(N)

for i in range(N):
    PPi = PPrange[i]
    offs = (EsatL / tauL).subs(v)
    # assume non-linear losses (qP(EP)) = 0:
    upperBound = (-offs + PPi * etaP / l).subs(v)
    # assume max. non-linear losses (qP(EP) = DR):
    lowerBound = (-offs + PPi * etaP / (l + DR)).subs(v)
    
    f = sym.lambdify(Pst, rootfunc.subs(PP, PPi))
    Pstrange[i] = brentq(f, lowerBound, upperBound)


    

In [4]:

    

%matplotlib inline
import matplotlib.pyplot as pl

# Initial slope efficiency (@ Pst = 0):
rho0 = (etaP * Toc / (l + DR)).subs(v)

fig, ax1 = pl.subplots()
pl.title('Output Power vs Pump Power')
pl.plot(PPrange, Toc.subs(v) * Pstrange, r'b-');
pl.plot([PPthreshold, PPrange[-1]], [0, rho0 * (PPrange[-1] - PPthreshold)],
        r'k--')
pl.xlabel('Pump Power (W)');
pl.ylabel('Output Power (W)', color='b');
ax2 = ax1.twinx()
ax2.plot(PPrange, Pstrange * (TR / EsatA).subs(v), r'r-')
ax2.set_ylabel('S-Parameter (1)', va='bottom', rotation=270, color='r');