{ SELF_FOCUS.PDE
This problem considers the self-focussing of a laser beam of Gaussian profile.
-- Submitted by John Trenholme
}
title
"2D GAUSSIAN BEAM PROFILE"
select
elevationgrid = 300
cubic { use cubic interpolation }
variables
real (threshold=0.1) { real (in-phase) part of field envelope }
imag (threshold=0.1) { imaginary (quadrature) part of field envelope }
definitions
radX = 2 { X "radius" of beam }
radY = 2 { Y "radius" of beam }
bMax = 2.25 { maximum B integral (= Time)}
zm = 5 { zoom-in factor for plots }
xHi = 7.17 * SQRT( radX * radY) { size of calculation domain... }
yHi = 7.17 * SQRT( radX * radY) { set for field = 0.001 at edge }
x45 = xHi * 0.7071 { point on boundary at 45 degrees }
y45 = yHi * 0.7071
tn =1e-30 { tiny number to force zero on plot scales }
power = PI * radX * radY * 2.73 ^ 2 / 8 { analytic integral }
inten = real * real + imag * imag { definition for later use }
time { "time" is really B integral }
0 to bMax by 0.03 * bMax
initial values
real = EXP( - ( x / ( radX * 2.73)) ^ 2 - ( y / ( radY * 2.73)) ^ 2)
imag = 0
equations { normalized, low-secular-phase nonlinear propagation }
real: DEL2( imag) + imag * ( inten - 1) = -DT( real)
imag: DEL2( real) + real * ( inten - 1) = DT( imag)
boundaries
region 1
start ( 0, 0) { bump is at center; only do one quadrant }
natural( real) = 0 { set slope to zero on boundary }
natural( imag) = 0 { if boundary value too big, move boundary out }
line to ( xHi, 0)
arc ( center = 0, 0) to ( 0, yHi)
line to ( 0, 0)
to close
monitors
for cycle = 1 { do this every cycle }
elevation( inten) from ( 0, 0) to ( xHi, 0) as "INTENSITY"
range( 0, tn)
contour( inten) as "INTENSITY" zoom ( 0, 0, xHi / zm, yHi / zm)
plots
for t = starttime { at the beginning only }
grid( x, y)
surface( inten) as "INTENSITY" range( 0, tn) viewpoint( 1000, 200, 40)
elevation( LOG10( inten)) from ( 0, 0) to ( xHi, 0) as "LOG10 INTENSITY"
for t = endtime { at the end only }
grid( x, y)
grid( x, y) zoom ( 0, 0, xHi / zm, yHi / zm)
surface( inten) as "INTENSITY" range( 0, tn) viewpoint( 1000, 200, 40)
zoom ( 0, 0, xHi / zm, yHi / zm)
elevation( LOG10( inten)) from ( 0, 0) to ( xHi, 0) as "LOG10 INTENSITY"
for t = starttime by ( endtime - starttime) / 5 to endtime { snapshots }
elevation( ARCTAN( imag / real) * 180 / PI) from ( 0, 0)
to ( xHi / zm, 0) as "PHASE (DEGREES)"
elevation( inten) from ( 0, 0) to ( xHi / zm, 0) as "INTENSITY"
range( 0, tn)
histories
history( inten) at ( 0, 0) ( 0.01 * xHi, 0) ( 0.03 * xHi, 0) ( 0.1 * xHi, 0)
( 0.3 * xHi, 0) ( xHi, 0) ( x45, y45) as "INTENSITY" print
history( real) at ( 0, 0) ( 0.01 * xHi, 0) ( 0.03 * xHi, 0) ( 0.1 * xHi, 0)
( 0.3 * xHi, 0) as "IN-PHASE FIELD"
history( imag) at ( 0, 0) ( 0.01 * xHi, 0) ( 0.03 * xHi, 0) ( 0.1 * xHi, 0)
( 0.3 * xHi, 0) as "QUADRATURE FIELD"
history( ARCTAN( imag / real) * 180 / PI) at ( 0, 0) ( 0.01 * xHi, 0)
( 0.03 * xHi, 0) ( 0.1 * xHi, 0) ( 0.3 * xHi, 0) as "PHASE (DEGREES)"
history( MIN( ( ABS( inten - 0.33)) ^ ( -0.75), 1)) at ( 0, 0)
range ( 0.045, 1) as "( INTENSITY - 0.33) ^ -0.75" { goes linearly to 0}
history( ABS( INTEGRAL( inten) / power - 1)) as "POWER CHANGE (EXACT = 0)"
end