\ slbench.fs
\
\ Benchmark version of sl.4th.
\
\ 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 "requires slbench 10 bench"
\	time gforth-fast slbench.fs -e "10 bench bye"

decimal

s" FLOATING-EXT" environment? [IF] ( flag) drop

\ 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-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 )
\	
\
\ ===============================
\ Handy definitions and constants
\ ===============================

: Sqr ( fx -- fx^2 ) fdup f* ;
: Mag2 ( fx fy -- fx^2+fy^2 ) Sqr fswap Sqr f+ ;

3 constant SVSIZE

: sv ( -- | defining word for a state vector )
\	create SVSIZE floats allot ;
	falign here SVSIZE floats allot constant ;

: sv! ( f1 f2 .. fn a -- | store a state vector at address a )
	SVSIZE 1- floats +
	\ SVSIZE 0 do dup >r f! r> 1 floats - loop drop ;
	SVSIZE 0 do dup f! 1 floats - loop drop ;

: sv@ ( a -- f1 f2 .. fn | fetch state vector from address a )
	\ SVSIZE 0 do dup >r f@ r> float+ loop drop ;
	SVSIZE 0 do dup f@ float+ loop drop ;

: sv+ ( a1 a2 -- f1 f2 .. fn | sum the state vectors at a1 and a2 and leave on stack )
	\ SVSIZE 0 do 2dup f@ rot f@ f+ 2swap float+ swap float+ swap loop 2drop ;
	SVSIZE 0 do 2dup f@ f@ f+ float+ swap float+ swap loop 2drop ;

: svc* ( f a -- f1 f2 .. fn | multiply the state vector at a by f )
	\ SVSIZE 0 do dup >r f@ fover f* fswap r> float+ loop drop fdrop ;
	SVSIZE 0 do dup f@ fover f* fswap float+ loop drop fdrop ;
  
-1e facos fconstant fpi

: intensity ( a -- fI | compute intensity of state vector )
	\ dup f@ rot float+ f@ Mag2 ;
	dup f@ float+ f@ Mag2 ;

: phase ( a -- fphase | compute phase in radians of state vector )
	\ dup f@ rot float+ f@ fswap fatan2 ;
	dup f@ float+ f@ fswap fatan2 ;

\ ================
\ Laser parameters
\ ================

fvariable t_p			\ photon lifetime (sec)
4.5e-12 t_p f!

fvariable t_s			\ carrier lifetime (sec)
700e-12 t_s f!

fvariable G_N			\ differential gain (cm^3/s)
2.6e-6 G_N f!

fvariable N_th			\ threshold carrier density (cm^-3)
1.5e18 N_th f!		

fvariable I_th			\ threshold current through laser (mA)
20e I_th f!

fvariable alpha			\ linewidth enhancement factor
5.0e alpha f!


\ ========================
\ Dimensionless parameters
\ ========================

fvariable T_ratio		\ T_ratio = t_s/t_p
fvariable PumpFactor 		\ P = PumpFactor*(I/I_th - 1)
fvariable P


\ =======================
\ Display all parameters
\ =======================

: 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! ;

: separator ( -- ) ." ===================================================" ;
: tab 9 emit ;
	
: params. ( -- | 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 
;  


init_params
\ params.

\ ======================================
\ Rates of dimensionless state variables
\ ======================================

\ Data in a state vector has the following order: Re{Y}, Im{Y}, Z

fvariable temp

: dY/ds ( a --  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} 	
;

fvariable temp2
		
: dZ/ds ( a  -- 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
	2e f* 1e 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/ ;


: dV/ds ( a -- 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 ;

\ =======================
\ ODE Solver
\ =======================

sv vstep

: stepV ( a ds ad -- 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 	
;


fvariable h
fvariable hh
fvariable h6

sv yt
sv dyt
sv dym


variable av
variable 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.
 
: rk4 ( a ds ad -- 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@ 6e 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!
	dyt dym sv+ dyt sv!
	adv @ dyt sv+ dyt sv!
	
	av @ h6 f@ dyt stepV
;


\ =============================
\ The injection current profile
\ =============================

fvariable fwhm		\ full-width at half-max for current pulse in ns
1e fwhm f!			( set to 1 ns )

fvariable pulse_amp	\ current pulse amplitude above d.c. level
20e pulse_amp f!		( set to 20 mA )

fvariable dc_current	\ d.c. current level
I_th f@ 10e f+ dc_current f!	( set to 10 mA above threshold )

fvariable peak_offset	\ offset in time for current peak
3e peak_offset f!		( set to 3 ns )

: GaussianPulse ( 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
	pulse_amp f@ f* dc_current f@ f+ ;		

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
\ =============================

: pump ( fc -- | compute and set the pump rate given the injection current )
	I_th f@ f/ 1e f-	\ compute P 
	PumpFactor f@ f* P f! ;

 
\ =============================
\ The rate equation solver
\ =============================

sv V			\ normalized state vector
sv Vdot			\ normalized time derivative of state vector

fvariable ds		\ dimensionless time step
0.1e ds f!		\ actual time step dt = t_p*ds


: main ( -- )

	init_params		\ compute all derived parameters
	2e fsqrt 0e 0e      	\ Re(E), Im(E), n 
	V sv!      		\ initialize the state vector

	\ Compute 20000 normalized time steps

	20000 0 do
	  i 0 d>f t_p f@ f* ds f@ f* 1e-9 f/ 	\ compute real-time in ns
	  \ fdup f. 2 spaces			\ output the time in ns
	  inj_current_function @ execute	\ compute the injection current
	  \ fdup f. 2 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
	    V intensity 
	    fdrop 
	    \ f. 2 spaces  \ output intensity
	    V phase fpi f/ 
	    fdrop 
	    \ f. 2 spaces  \ output normalized phase: 1.0 = pi radians
	    V 2 floats + f@ 
	    fdrop 
	    \ f. cr \ output normalized carrier density 
	  then
	loop
;

\ printing version for debugging
: sl. ( #timesteps -- )

	init_params		\ compute all derived parameters
	2e fsqrt 0e 0e      	\ Re(E), Im(E), n 
	V sv!      		\ initialize the state vector

	\ Compute 20000 normalized time steps

	\ 20000 0 do
        cr
	( #timesteps) 0 do
	  i 0 d>f t_p f@ f* ds f@ f* 1e-9 f/ 	\ compute real-time in ns
	  fdup f. 2 spaces			\ output the time in ns
	  inj_current_function @ execute	\ compute the injection current
	  fdup f. 2 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
	    V intensity 
	    \ fdrop 
	    f. 2 spaces  \ output intensity
	    V phase fpi f/ 
	    \ fdrop 
	    f. 2 spaces  \ output normalized phase: 1.0 = pi radians
	    V 2 floats + f@ 
	    \ fdrop 
	    f. cr \ output normalized carrier density 
	  then
	loop
;

\ =============================
\ Run benchmark
\ =============================

: bench  ( u -- )  ( u) 0 DO main LOOP ;

[ELSE] .( ***Floating point not available.) [THEN]