\ slbench.hf \ \ ^Forth/pfe version of slbench.fs (benchmark version of sl.4th) \ for automatic translation to C. \ \ Copyright (c) 2000--2002 Krishna Myneni \ \ This library is free software; you can redistribute it and/or \ modify it under the terms of the GNU Lesser General Public \ License as published by the Free Software Foundation; either \ version 2.1 of the License, or (at your option) any later \ version. \ This library is distributed in the hope that it will be useful, \ but WITHOUT ANY WARRANTY; without even the implied warranty of \ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the \ GNU Lesser General Public License for more details. \ \ Usage: \ \ time pfe -q -y -e "loadm floating loadm slbench 10 bench" \ s" FLOATING-EXT" environment? \ 0= [IF] .( ***Floating point not available.) \ [ELSE] ( flag) drop [THEN] \ decimal \ Solve the semiconductor laser rate equations, given a current pulse \ profile. Output the time, intensity, phase, and current density. \ \ Based on the program ned10.c by S.D. Pethel \ \ Krishna Myneni, 1-26-2000 \ Further modifications for portability by David N. Williams, 2002. \ \ \ Revisions: \ \ 2002-12-6 modified SV to ensure portability of float \ alignment. DNW \ \ 2002-11-19 added printing version of main loop -- execute \ <#timesteps> sl. DNW \ \ 2002-11-07 modified for ^Forth. DNW \ \ 2002-10-27 changed all instances of dfloat to float for \ ANS Forth portability. Removed explicit fp \ number size dep. KM \ \ 2002-10-24 fixed problem with the main loop; prev. was not \ computing Vdot on every loop iteration. Also \ changed current pulse pos to 3 ns. KM \ \ 2002-10-21 fixed time scale problem in sl after problem \ was pointed out by Marcel Hendrix. KM \ \ 5-12-2001 this version modified to work on Forth systems \ which use a separate floating point stack. Output \ to console has been commented out in "main" -- KM \ \ ------------------------------------------------------------------- \ \ The normalized laser rate equations are given by \ (cf D.W. Sukow, PhD Thesis: Experimental Control of \ Instabilities and Chaos in Fast Dynamical Systems, 1997): \ \ dY/ds = (1 + i*alpha)ZY \ dZ/ds = 1/T (P - Z - (1 + 2Z)|Y|^2) \ \ obtained by the transformation \ \ Y = (t_s*G_N/2)^0.5 * E ( note E is complex ) \ \ Z = (t_p*G_N/2)(N - N_th) \ \ s = t/t_p ( time in units of photon lifetime ) \ \ and the parameters are given by: \ \ P = (t_p*G_N*N_th/2)(I/I_th - 1) \ T = t_s/t_p \ \ The basic quantities are: \ \ E ( complex electric field in photons^0.5/cm^3 ) \ N ( carrier density in cm^-3 ) \ N_th ( carrier density at threshold for lasing ) \ I_th ( threshold current in mA ) \ t_p ( photon lifetime in sec ) \ t_s ( carrier lifetime in sec ) \ alpha ( linewidth enhancement factor -- no dimensions ) \ G_N ( differential gain at threshold in cm^3/s ) \ \ BEGIN-EXTERNAL {" /* Extra includes. */ /* #include */ #include "} REQUIRES hfkernel.hf REQUIRES hffloat.hf END-EXTERNAL MAKE-INDEX slbench "KM's laser rate benchmark" \ =============================== \ Handy definitions and constants \ =============================== \ target magic 0e fconstant: fzero "f_zero" 1e fconstant: fone "f_one" 2e fconstant: ftwo "f_two" 6e fconstant: fsix "f_six" 0.5e fconstant: fhalf "f_half" def: Sqr "sqr" (f: fx -- fx^2 ) fdup f* ;def def: Mag2 "mag2" (f: fx fy -- fx^2+fy^2 ) Sqr fswap Sqr f+ ;def host[ 3 constant [SVSIZE] ]host [SVSIZE] constant: SVSIZE "SVSIZE" host[ : sv ( -- | defining word for a state vector ) \ falign here SVSIZE floats allot constant ; [SVSIZE] fcreate-allot: ; ]host def: sv! "sv_store" (f: f1 f2 .. fn s: a -- | store a state vector at address a ) SVSIZE 1- floats + \ SVSIZE 0 do dup f! 1 floats - loop drop ; SVSIZE zero do dup f! one floats - loop drop ;def def: sv@ "sv_fetch" ( a -- f: f1 f2 .. fn | fetch state vector from address a ) \ SVSIZE 0 do dup f@ float+ loop drop ; SVSIZE zero do dup f@ float+ loop drop ;def def: sv+ "sv_plus" ( a1 a2 -- f: f1 f2 .. fn | sum the state vectors at a1 and a2 and leave on stack ) \ SVSIZE 0 do 2dup f@ f@ f+ float+ swap float+ swap loop 2drop ; SVSIZE zero do 2dup f@ f@ f+ float+ swap float+ swap loop 2drop ;def def: svc* "svc_star" ( a f: f-- f1 f2 .. fn | multiply the state vector at a by f ) \ SVSIZE 0 do dup f@ fover f* fswap float+ loop drop fdrop ; SVSIZE zero do dup f@ fover f* fswap float+ loop drop fdrop ;def host[ -1e facos ]host fconstant: fpi "fpi" def: intensity "intensity" ( a -- f: fI | compute intensity of state vector ) \ dup f@ rot float+ f@ Mag2 ; dup f@ float+ f@ Mag2 ;def def: phase "phase" ( a -- fphase | compute phase in radians of state vector ) \ dup f@ rot float+ f@ fswap fatan2 ; dup f@ float+ f@ fswap fatan2 ;def \ ================ \ Laser parameters \ ================ 4.5e-12 init-fvariable: t_p "t_p" \ photon lifetime (sec) \ 4.5e-12 t_p f! 700e-12 init-fvariable: t_s "t_s" \ carrier lifetime (sec) \ 700e-12 t_s f! 2.6e-6 init-fvariable: G_N "G_N" \ differential gain (cm^3/s) \ 2.6e-6 G_N f! 1.5e18 init-fvariable: N_th "N_th" \ threshold carrier density (cm^-3) \ 1.5e18 N_th f! host[ 20e fconstant I_th_0 ]host I_th_0 init-fvariable: I_th "I_th" \ threshold current through laser (mA) \ 20e I_th f! 5.0e init-fvariable: alpha "alpha" \ linewidth enhancement factor \ 5.0e alpha f! \ ======================== \ Dimensionless parameters \ ======================== fvariable: T_ratio "T_ratio" \ T_ratio = t_s/t_p fvariable: PumpFactor "PumpFactor" \ P = PumpFactor*(I/I_th - 1) fvariable: P "P" \ ======================= \ Display all parameters \ ======================= def: init_params "init_params" ( -- | compute the normalized parameters ) t_s f@ t_p f@ f/ T_ratio f! \ t_p f@ G_N f@ f* N_th f@ f* 2e f/ PumpFactor f! ; t_p f@ G_N f@ f* N_th f@ f* ftwo f/ PumpFactor f! ;def def: separator "separator" ( -- ) ." ===================================================" ;def def: tab "tab" nine emit ;def def: params. "params_dot" ( -- | display all of the parameters ) cr separator cr ." Symbol" tab ." Parameter " tab ." Value" cr separator cr cr ." t_p " tab ." Photon lifetime (s): " tab t_p f@ f. cr ." t_s " tab ." Carrier lifetime (s): " tab t_s f@ f. cr ." G_N " tab ." Differential gain (cm^3/s): " tab G_N f@ f. cr ." N_th " tab ." Thr. carrier density (cm^-3): " tab N_th f@ f. cr ." I_th " tab ." Thr. current (mA): " tab I_th f@ f. cr ." alpha" tab ." Linewidth enhancement factor: " tab alpha f@ f. cr separator cr ." Derived Dimensionless Parameters " cr separator cr cr ." t_s/t_p ratio: " tab T_ratio f@ f. cr ." Pump factor: " tab PumpFactor f@ f. cr separator cr ;def def: init_params. "init_params_dot" ( -- ) init_params params. ;def \ ====================================== \ Rates of dimensionless state variables \ ====================================== \ Data in a state vector has the following order: Re{Y}, Im{Y}, Z fvariable: temp "temp" def: dY/ds "dY_by_ds" ( a -- f: Re{dY/ds} Im{dY/ds} | compute dY/ds for state vector 'a' ) \ Re{dY/ds} = Z(Re{Y} - alpha*Im{Y}) \ \ Im{dY/ds} = Z(alpha*Re{Y} + Im{Y}) sv@ \ fdup 2>r 2>r temp f! \ store z fover fover \ alpha f@ f* f- 2r> f* alpha f@ f* f- temp f@ f* \ Re{dY/ds} frot alpha f@ f* \ frot f+ 2r> f* frot f+ temp f@ f* \ Im{dY/ds} ;def fvariable: temp2 "temp2" def: dZ/ds "dZ_by_ds" ( a -- f: dZ/ds | compute rate of change of Z ) \ dZ/ds = (P - Z - (1 + 2Z)|Y|^2)/T \ \ P = PumpFactor*(I/I_th - 1); P should be computed prior to executing dZ/ds sv@ \ fdup 2>r fdup temp f! \ 2e f* 1e f+ 2>r ftwo f* fone f+ temp2 f! \ Mag2 2r> f* Mag2 temp2 f@ f* \ P f@ 2r> f- fswap f- P f@ temp f@ f- fswap f- T_ratio f@ f/ ;def def: dV/ds "dV_by_ds" ( a -- f: Re{dY/ds} Im{dY/ds} dZ/ds | derivative of the state vector 'a' ) \ dup dZ/ds 2>r dY/ds 2r> ; dup dZ/ds dY/ds frot ;def \ ======================= \ ODE Solver \ ======================= sv vstep "vstep" def: stepV "stepV" ( a ad f:ds -- Re{Y}' Im{Y}' Z' | compute stepped vector ) \ Compute V' = V + ds*(dV/ds) \ a is the address of the state vector V \ ds is the step size (floating point) \ ad is the address of a vector containing dV/ds svc* vstep sv! \ multiply derivatives by ds to obtain the step vstep sv+ \ add step to initial vector ;def fvariable: h "h" fvariable: hh "hh" fvariable: h6 "h6" sv yt "yt" sv dyt "dyt" sv dym "dym" variable: av "av" variable: adv "adv" \ rk4 is an *approximate* fourth-order Runge-Kutta ODE stepper, \ specific to this problem. This stepper is only valid when \ the derivatives dV/ds do not change appreciably over a \ single step. It ignores the explicit time dependence \ of the derivative at internal points in the time step. \ This assumption is valid for the case here because the \ current does not vary much over one time step of ds*t_p. def: rk4 "rk4" ( a ad f: ds -- f1 ... fn | final state vector is returned on the stack ) \ a is the address of the state vector \ ds is the step size \ ad is the address of a vector containing the initial derivatives adv ! h f! av ! \ h f@ 0.5e f* hh f! \ hh = h/2; h f@ fhalf f* hh f! \ hh = h/2; \ h f@ 6e f/ h6 f! \ h6 = h/6; h f@ fsix f/ h6 f! \ h6 = h/6; av @ hh f@ adv @ stepV \ first step to midpoint yt sv! yt dV/ds dyt sv! \ second step av @ hh f@ dyt stepV yt sv! yt dV/ds dym sv! \ third step av @ h f@ dym stepV yt sv! dyt dym sv+ dym sv! yt dV/ds dyt sv! \ fourth step \ 2e dym svc* dym sv! ftwo dym svc* dym sv! dyt dym sv+ dyt sv! adv @ dyt sv+ dyt sv! av @ h6 f@ dyt stepV ;def \ ============================= \ The injection current profile \ ============================= 1e init-fvariable: fwhm "fwhm" \ full-width at half-max for current pulse in ns \ 1e fwhm f! ( set to 1 ns ) 20e init-fvariable: pulse_amp "pulse_amp" \ current pulse amplitude above d.c. level \ 20e pulse_amp f! ( set to 20 mA ) host[ I_th_0 10e f+ ]host init-fvariable: dc_current "d_c_current" \ d.c. current level \ I_th f@ 10e f+ dc_current f! ( set to 10 mA above threshold ) 3e init-fvariable: peak_offset "peak_offset" \ offset in time for current peak \ 3e peak_offset f! ( set to 3 ns ) -2.77066e fconstant: Pulse_magic "Pulse_magic" def: GaussianPulse "GaussianPulse" (f: ft -- fc | compute current at real time ft ) \ ft is in nano-seconds peak_offset f@ f- fwhm f@ f/ \ Sqr -2.77066e f* fexp Sqr Pulse_magic f* fexp pulse_amp f@ f* dc_current f@ f+ ;def \ Have to think about how to handle EXECUTE in ^Forth... \ variable inj_current_function \ address of the injection current generating function \ ' GaussianPulse inj_current_function ! \ You may use your own injection current word by typing \ \ ' yourword inj_current_function ! \ \ The word should have the stack diagram ( ft -- fc ) where ft \ is the time in nanoseconds, and fc is the current in mA. \ ============================= \ Pump rate \ ============================= def: pump "pump" (f: fc -- | compute and set the pump rate given the injection current ) \ I_th f@ f/ 1e f- \ compute P I_th f@ f/ fone f- \ compute P PumpFactor f@ f* P f! ;def \ ============================= \ The rate equation solver \ ============================= sv V "V" \ normalized state vector sv Vdot "Vdot" \ normalized time derivative of state vector 0.1e init-fvariable: ds "ds" \ dimensionless time step \ 0.1e ds f! \ actual time step dt = t_p*ds \ target magic 20000 constant: #time-steps "sharp_time_steps" \ 100 constant: #time-steps "sharp_time_steps" 1e-9 fconstant: ns "ns" def: main "main" ( -- ) init_params \ compute all derived parameters \ 2e fsqrt 0e 0e \ Re(E), Im(E), n ftwo fsqrt fzero fzero \ Re(E), Im(E), n V sv! \ initialize the state vector \ Compute 20000 normalized time steps \ 20000 0 do #time-steps zero do \ i 0 d>f t_p f@ f* ds f@ f* 1e-9 f/ \ compute real-time in ns i zero d>f t_p f@ f* ds f@ f* ns f/ \ compute real-time in ns \ \ fdup f. 2 spaces \ output the time in ns \ fdup f. two spaces \ output the time in ns \ inj_current_function @ execute \ compute the injection current GaussianPulse \ \ fdup f. 2 spaces \ output the current in mA \ fdup f. two spaces \ output the current in mA pump V dV/ds Vdot sv! \ evaluate derivatives at current time V ds f@ Vdot rk4 \ step the state vector V sv! \ i 1 mod 0= if i one mod 0= if V intensity fdrop \ \ f. 2 spaces \ output intensity \ f. two spaces \ output intensity V phase fpi f/ fdrop \ \ f. 2 spaces \ output normalized phase: 1.0 = pi radians \ f. two spaces \ output normalized phase: 1.0 = pi radians \ V 2 floats + f@ V two floats + f@ fdrop \ f. cr \ output normalized carrier density then loop ;def \ Printing version for debugging. def: sl. "s_l_dot" ( #timesteps -- ) init_params \ compute all derived parameters ftwo fsqrt fzero fzero \ Re(E), Im(E), n V sv! \ initialize the state vector \ Compute 20000 normalized time steps cr ( #time-steps) zero do i zero d>f t_p f@ f* ds f@ f* ns f/ \ compute real-time in ns fdup f. two spaces \ output the time in ns GaussianPulse \ compute the injection current fdup f. two spaces \ output the current in mA pump V dV/ds Vdot sv! \ evaluate derivatives at current time V ds f@ Vdot rk4 \ step the state vector V sv! i one mod 0= if V intensity \ fdrop f. two spaces \ output intensity V phase fpi f/ \ fdrop f. two spaces \ output normalized phase: 1.0 = pi radians V two floats + f@ \ fdrop f. cr \ output normalized carrier density then loop ;def \ ============================= \ Run benchmark \ ============================= def: bench "bench" ( u -- ) ( u) zero DO main LOOP ;def