# Rcjp's Weblog

## April 4, 1993

### Final Year BSc Project

Filed under: physics — rcjp @ 7:32 pm

I did some initial work on the effects of soliton propagation in birefringent fibres around this time writing a program I later cleaned up a bit for my PhD.

There were a few dozen (literally) versions of this lying around on my HDD, back then you had to write a new program and compile it for each new set of parameters you were investigating – which meant leaving lots of versions of code around.

Its one of the things that drove me to buy Mathematica 2.0 (I think it was). I developed a system for splicing the output from mathematica into a fortran program and running that – but it wasn’t general enough to be useful. The mathematica syntax, which I struggled using, led me to Scheme and finally to Common Lisp – I wish I could have used it back then!

Essentially I wrote a program to solve:
${\partial u\over\partial\xi}=-\delta{\partial u\over\partial\tau}+{i\over 2}{\partial^2 u\over\partial\tau^2}+i\left[|u|^2+{2\over 3}|v|^2\right] u$
${\partial v\over\partial\xi}=+\delta{\partial v\over\partial\tau}+{i\over 2}{\partial^2 v\over\partial\tau^2}+i\left[|v|^2+{2\over 3}|u|^2\right] v$

where u and v are different modes of the fibre. Fortran code…

    CCCC===================================================================
C  LAUNCHES A SOLITON OF VARYING AMPLITUDE 0.4, 0.8 AT THETA = 45
CCCC===================================================================

PROGRAM SPLSTX
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION U01(1024),UZ1(1024),UZH1(1024),U01S(1024),
*          U02(1024),UZ2(1024),UZH2(1024),U02S(1024)
DOUBLE PRECISION N
PARAMETER (THDS = 2.0D0/3.0D0)

OPEN (10,FILE = 'RES1AA')
OPEN (11,FILE = 'RES1AB')
OPEN (12,FILE = 'RES1BA')
OPEN (13,FILE = 'RES1BB')
OPEN (14,FILE = 'RES1CA')
OPEN (15,FILE = 'RES1CB')
OPEN (16,FILE = 'RES1DA')
OPEN (17,FILE = 'RES1DB')
OPEN (18,FILE = 'RES1U')
OPEN (19,FILE = 'RES1V')

OPEN (20,FILE = 'RES2AA')
OPEN (21,FILE = 'RES2AB')
OPEN (22,FILE = 'RES2BA')
OPEN (23,FILE = 'RES2BB')
OPEN (24,FILE = 'RES2CA')
OPEN (25,FILE = 'RES2CB')
OPEN (26,FILE = 'RES2DA')
OPEN (27,FILE = 'RES2DB')
OPEN (28,FILE = 'RES2U')
OPEN (29,FILE = 'RES2V')

DO 2000 MNUMBER = 1,2
N = MNUMBER*0.4D0

CCCC------------------------------------------------------------------
CCCC N      - ORDER OF SOLITON
CCCC NUMZ   - NUMBER OF SLICES
CCCC NOWRIT - NUMBER OF ITERATIONS WHERE NO RESULT IS WRITTEN OUT
CCCC NN     - NUMBER OF POINTS ACROSS PULSE
CCCC RT     - WIDTH OF PULSE
CCCC H      - STEPSIZE BETWEEN SLICES
CCCC THETA  - ANGLE AT WHICH SOLITON IS LAUNCHED
CCCC VWIND  - VIEW WINDOW OF PULSE
CCCC------------------------------------------------------------------

NUMZ    = 21
NOWRIT  = 250
NN      = 512
THETA   = 45.0D0
RT      = 220.0D0
DELTA   = 0.5D0
VWIND   = 60.0d0

PI    = 4.0D0*DATAN(1.0D0)
RTO2  = RT/2.0D0
RTONN = RT/NN
NNO2  = NN/2.0D0
H     = 30.0D0*PI/(NUMZ*NOWRIT)
COEFF = DBLE(NN)
SLOW  = DCOS(THETA*PI/180.0D0)
FAST  = DSIN(THETA*PI/180.0D0)

CCCC CLEAR ARRAYS

CALL ZERO_R(U01,NN*2)
CALL ZERO_R(U01S,NN*2)
CALL ZERO_R(UZ1,NN*2)
CALL ZERO_R(UZH1,NN*2)

CALL ZERO_R(U02,NN*2)
CALL ZERO_R(U02S,NN*2)
CALL ZERO_R(UZ2,NN*2)
CALL ZERO_R(UZH2,NN*2)

CCCC GENERATE FIRST PULSE

DO I=1,2*NN,2
TIME    = -RTO2+(I/2-1)*RTONN
U01(I)  = N*SLOW/DCOSH(TIME)
U01S(I) = U01(I)*U01(I)
ENDDO

CCCC GENERATE SECOND PULSE

DO I=1,2*NN,2
TIME    = -RTO2+(I/2-1)*RTONN
U02(I)  = N*FAST/DCOSH(TIME)
U02S(I) = U02(I)*U02(I)
ENDDO

CCCC COPY INITIAL PULSE TO UZ() WHICH IS USED TO
CCCC CALCULATED THE PULSE SHAPE AT Z+H/2
CCCC AND UZH

DO 10 I=1,2*NN
UZ1(I)  = U01(I)
UZH1(I) = U01(I)
UZ2(I)  = U02(I)
UZH2(I) = U02(I)
10   CONTINUE

CCCC --------START OF MAIN LOOP - PROPAGATES PULSE---------

DO 1000 IYCOUNT=1,NUMZ

CCCC WRITE PULSE TO FILE
if (iycount.eq.1) then
DO I=1,2*NN,2
TIME = -RTO2+(I/2-1)*RTONN
if (abs(time).lt.vwind) then
write (MNUMBER*10,*) time, SQRT(U01S(I))
write (MNUMBER*10+1,*) time, SQRT(U02S(I))
endif
ENDDO
endif

if (iycount.eq.7) then
DO I=1,2*NN,2
TIME = -RTO2+(I/2-1)*RTONN
if (abs(time).lt.vwind) then
write (MNUMBER*10+2,*) time, SQRT(U01S(I))
write (MNUMBER*10+3,*) time, SQRT(U02S(I))
endif
ENDDO
endif

if (iycount.eq.14) then
DO I=1,2*NN,2
TIME = -RTO2+(I/2-1)*RTONN
if (abs(time).lt.vwind) then
write (MNUMBER*10+4,*) time, SQRT(U01S(I))
write (MNUMBER*10+5,*) time, SQRT(U02S(I))
endif
ENDDO
endif

if (iycount.eq.21) then
DO I=1,2*NN,2
TIME = -RTO2+(I/2-1)*RTONN
if (abs(time).lt.vwind) then
write (MNUMBER*10+6,*) time, SQRT(U01S(I))
write (MNUMBER*10+7,*) time, SQRT(U02S(I))
endif
ENDDO
endif

DO I=1,2*NN,2
TIME = -RTO2+(I/2-1)*RTONN
if (abs(time).lt.vwind) then
write (MNUMBER*10+8,*) time, SQRT(U01S(I))
write (MNUMBER*10+9,*) time, SQRT(U02S(I))
endif
ENDDO

PRINT '(20X,I4,3X,A,3X,I4)',IYCOUNT,'OF',NUMZ
PRINT *,' M NUMBER IS ',MNUMBER

CCCC DUMMY LOOP PROPAGATES PULSE WITHOUT WRITING OUT RESULTS

DO 500 IDUMMY=1,NOWRIT

CCCC CALCULATE DISPERSIVE EFFECTS BETWEEN Z -> Z + H/2 FOR 1ST PULSE

ISIGN = 1
CALL FOUR1(UZ1,NN,ISIGN)

DO I=1,NN+1,2
W1       = 2.0D0*PI/RT*((I-1)/2)
THETA    = -0.5D0*H*(DELTA*W1+W1*W1*0.5D0)
TEMP     = UZ1(I)
C        = DCOS(THETA)
S        = DSIN(THETA)
UZ1(I)   = UZ1(I)*C-UZ1(I+1)*S
UZ1(I+1) = TEMP*S+UZ1(I+1)*C
ENDDO

DO I=NN+3,NN*2-1,2
W1       = 2.0D0*PI/RT*(-NN+(I-1)/2)
THETA    = -0.5D0*H*(DELTA*W1+W1*W1*0.5D0)
TEMP     = UZ1(I)
C        = DCOS(THETA)
S        = DSIN(THETA)
UZ1(I)   = UZ1(I)*C-UZ1(I+1)*S
UZ1(I+1) = TEMP*S+UZ1(I+1)*C
ENDDO

ISIGN = -1
CALL FOUR1(UZ1,NN,ISIGN)

CCCC CALCULATE DISPERSIVE EFFECTS BETWEEN Z -> Z + H/2 FOR 2ND PULSE

ISIGN = 1
CALL FOUR1(UZ2,NN,ISIGN)

DO I=1,NN+1,2
W2       = 2.0D0*PI/RT*((I-1)/2)
THETA    = 0.5D0*H*(DELTA*W2-W2*W2*0.5D0)
C        = DCOS(THETA)
S        = DSIN(THETA)
TEMP     = UZ2(I)
UZ2(I)   = UZ2(I)*C-UZ2(I+1)*S
UZ2(I+1) = TEMP*S+UZ2(I+1)*C
ENDDO

DO I=NN+3,NN*2-1,2
W2       = 2.0D0*PI/RT*(-NN+(I-1)/2)
THETA    = 0.5D0*H*(DELTA*W2-W2*W2*0.5D0)
TEMP     = UZ2(I)
C        = DCOS(THETA)
S        = DSIN(THETA)
UZ2(I)   = UZ2(I)*C-UZ2(I+1)*S
UZ2(I+1) = TEMP*S+UZ2(I+1)*C
ENDDO

ISIGN = -1
CALL FOUR1(UZ2,NN,ISIGN)

CCCC MULTIPLY BOTH RESULTS BY 1/NN

DO 60 I=1,2*NN
UZ1(I)  = UZ1(I)/COEFF
UZ2(I)  = UZ2(I)/COEFF
60        CONTINUE

CCCC ---------CALCULATE NONLINEAR EFFECTS------------

CCCC J-LOOP USED TO FIND THE PULSE SHAPE AT Z+H
CCCC STARTING WITH THE APPROXIMATION THAT IT IS THE PULSE IN Z

DO 200 J=1,3
DO 70 I=1,2*NN,2
CCCC FIRST PULSE
THETA =  0.5D0*H*
*                   ( U01S(I) + THDS*U02S(I) +
*                     UZH1(I)*UZH1(I) + UZH1(I+1)*UZH1(I+1) +
*               THDS*(UZH2(I)*UZH2(I) + UZH2(I+1)*UZH2(I+1)) )
C1  = DCOS(THETA)
S1  = DSIN(THETA)
CCCC SECOND PULSE
THETA = 0.5D0*H*
*                   ( U02S(I) + THDS*U01S(I) +
*                     UZH2(I)*UZH2(I) + UZH2(I+1)*UZH2(I+1) +
*               THDS*(UZH1(I)*UZH1(I) + UZH1(I+1)*UZH1(I+1)) )
C2  = DCOS(THETA)
S2  = DSIN(THETA)

CCCC MULTIPLY WITH DISPERSION ONLY TERM UZ()

UZH1(I)   = UZ1(I)*C1 - UZ1(I+1)*S1
UZH1(I+1) = UZ1(I)*S1 + UZ1(I+1)*C1
UZH2(I)   = UZ2(I)*C2 - UZ2(I+1)*S2
UZH2(I+1) = UZ2(I)*S2 + UZ2(I+1)*C2
70          CONTINUE

CCCC CALCULATE DISPERSIVE EFFECTS BETWEEN Z + H/2 -> Z + H 1ST PULSE

ISIGN = 1
CALL FOUR1(UZH1,NN,ISIGN)

DO I=1,NN+1,2
W1        = 2.0D0*PI/RT*((I-1)/2)
THETA    = -0.5D0*H*(DELTA*W1+W1*W1*0.5D0)
TEMP      = UZH1(I)
C         = DCOS(THETA)
S         = DSIN(THETA)
UZH1(I)   = UZH1(I)*C-UZH1(I+1)*S
UZH1(I+1) = TEMP*S+UZH1(I+1)*C
ENDDO

DO I=NN+3,NN*2-1,2
W1        = 2.0D0*PI/RT*(-NN+(I-1)/2)
THETA    = -0.5D0*H*(DELTA*W1+W1*W1*0.5D0)
TEMP      = UZH1(I)
C         = DCOS(THETA)
S         = DSIN(THETA)
UZH1(I)   = UZH1(I)*C-UZH1(I+1)*S
UZH1(I+1) = TEMP*S+UZH1(I+1)*C
ENDDO

ISIGN = -1
CALL FOUR1(UZH1,NN,ISIGN)

CCCC CALCULATE DISPERSIVE EFFECTS BETWEEN Z + H/2 -> Z + H 2ND PULSE

ISIGN = 1
CALL FOUR1(UZH2,NN,ISIGN)

DO I=1,NN+1,2
W2        = 2.0D0*PI/RT*((I-1)/2)
THETA     = 0.5D0*H*(DELTA*W2-W2*W2*0.5D0)
TEMP      = UZH2(I)
C         = DCOS(THETA)
S         = DSIN(THETA)
UZH2(I)   = UZH2(I)*C-UZH2(I+1)*S
UZH2(I+1) = TEMP*S+UZH2(I+1)*C
ENDDO

DO I=NN+3,NN*2-1,2
W2        = 2.0D0*PI/RT*(-NN+(I-1)/2)
THETA     = 0.5D0*H*(DELTA*W2-W2*W2*0.5D0)
TEMP      = UZH2(I)
C         = DCOS(THETA)
S         = DSIN(THETA)
UZH2(I)   = UZH2(I)*C-UZH2(I+1)*S
UZH2(I+1) = TEMP*S+UZH2(I+1)*C
ENDDO

ISIGN = -1
CALL FOUR1(UZH2,NN,ISIGN)

CCCC MULTIPLY BY 1/NN BOTH PULSES

DO 120 I=1,2*NN
UZH1(I) = UZH1(I)/COEFF
UZH2(I) = UZH2(I)/COEFF
120          CONTINUE
200       CONTINUE

CCCC RE-INITIALISE U0(), UZ() AND UZH() TO BE THE NEW PULSE SHAPE

DO 220 I=1,2*NN
UZ1(I)  = UZH1(I)
U01(I)  = UZ1(I)

UZ2(I)  = UZH2(I)
U02(I)  = UZ2(I)
220       CONTINUE
DO I=1,2*NN,2
U01S(I) = U01(I)*U01(I)+U01(I+1)*U01(I+1)
U02S(I) = U02(I)*U02(I)+U02(I+1)*U02(I+1)
ENDDO

500    CONTINUE
1000 CONTINUE
2000 CONTINUE
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(13)
CLOSE(14)
CLOSE(15)
CLOSE(16)
CLOSE(17)
CLOSE(18)
CLOSE(19)
CLOSE(20)
CLOSE(21)
CLOSE(22)
CLOSE(23)
CLOSE(24)
CLOSE(25)
CLOSE(26)
CLOSE(27)
CLOSE(28)
CLOSE(29)
END

SUBROUTINE ZERO_R(U,NN)
CCCC==================================================================
C    CLEAR THE ARRAYS
CCCC==================================================================
DOUBLE PRECISION U(*)
DO 10 I=1,NN
U(I) = 0.0D0
10 CONTINUE
END

SUBROUTINE FOUR1(DATAM,NN,ISIGN)
CCCC==================================================================
C    FOURIER TRANSFORMS DATAM ARRAY (ISIGN=1) INVERSE (ISIGN=-1)
C    DATAM ARRAY CONTAINS ALTERNATE REAL AND IMAGINARY PAIRS
C    NN MUST BE AN INTEGER POWER OF 2
CCCC==================================================================
DOUBLE PRECISION WR,WI,WPR,WPI,WTEMP,THETA
DOUBLE PRECISION DATAM(2*NN)
DOUBLE PRECISION TEMPR,TEMPI
INTEGER ISIGN,NN,I,ISTEP,J,M,MMAX,N
N = 2*NN
J = 1
DO 11 I=1,N,2
IF (J.GT.I) THEN
TEMPR      = DATAM(J)
TEMPI      = DATAM(J+1)
DATAM(J)   = DATAM(I)
DATAM(J+1) = DATAM(I+1)
DATAM(I)   = TEMPR
DATAM(I+1) = TEMPI
ENDIF
M = N/2
1    IF ((M.GE.2).AND.(J.GT.M)) THEN
J = J-M
M = M/2
GOTO 1
ENDIF
J = J+M
11 CONTINUE
MMAX = 2
2 IF (N.GT.MMAX) THEN
ISTEP = 2*MMAX
THETA = 6.28318530717959D0/(ISIGN*MMAX)
WPR   = -2.D0*DSIN(0.5D0*THETA)**2
WPI   = DSIN(THETA)
WR    = 1.D0
WI    = 0.D0
DO 13 M=1,MMAX,2
DO 12 I=M,N,ISTEP
J          = I+MMAX
TEMPR      = WR*DATAM(J)-WI*DATAM(J+1)
TEMPI      = WR*DATAM(J+1)+WI*DATAM(J)
DATAM(J)   = DATAM(I)-TEMPR
DATAM(J+1) = DATAM(I+1)-TEMPI
DATAM(I)   = DATAM(I)+TEMPR
DATAM(I+1) = DATAM(I+1)+TEMPI
12      CONTINUE
WTEMP = WR
WR    = WR*WPR-WI*WPI+WR
WI    = WI*WPR+WTEMP*WPI+WI
13   CONTINUE
MMAX = ISTEP
GOTO 2
ENDIF
RETURN
END


and I used GLE to get some quite nice graphs

    size  21 29
amove 3 19.3
begin graph
size 8 9
nobox
title "N = 0.4" hei 0.6 font texcmr
ytitle "Amplitude" hei 0.5
xlabels off
!        yaxis min 0.0
data res1aa
data res1ab
d1 lstyle 1 smooth
d2 lstyle 2 smooth
end graph
rmove 0 -6.5
begin graph
size 8 9
nobox
xlabels off
ytitle "Amplitude" hei 0.5
!        yaxis min 0.0
data res1ba
data res1bb
d1 lstyle 1 smooth
d2 lstyle 2 smooth
end graph
rmove  0 -6.5
begin graph
size 8 9
nobox
xlabels off
ytitle "Amplitude" hei 0.5
!        yaxis min 0.0
data res1ca
data res1cb
d1 lstyle 1 smooth
d2 lstyle 2 smooth
end graph
rmove  0 -6.5
begin graph
size 8 9
nobox
xtitle "Time \tau" hei 0.5
ytitle "Amplitude" hei 0.5
!        yaxis min 0.0
data res1da
data res1db
d1 lstyle 1 smooth
d2 lstyle 2 smooth
end graph
! second graph
rmove 8 19.5
begin graph
size 8 9
nobox
title "N = 0.8" hei 0.6 font texcmr
xlabels off
y2title "\xi = 0\pi" hei 0.6 font texcmr rotate
data res2aa
data res2ab
d1 lstyle 1 smooth
d2 lstyle 2 smooth
end graph
rmove 0 -6.5
begin graph
size 8 9
nobox
xlabels off
y2title "\xi = 10\pi" hei 0.6 font texcmr rotate
data res2ba
data res2bb
d1 lstyle 1 smooth
d2 lstyle 2 smooth
end graph
rmove  0 -6.5
begin graph
size 8 9
nobox
xlabels off
y2title "\xi = 20\pi" hei 0.6 font texcmr rotate
data res2ca
data res2cb
d1 lstyle 1 smooth
d2 lstyle 2 smooth
end graph
rmove  0 -6.5
begin graph
size 8 9
nobox
xtitle "Time \tau" hei 0.5
y2title "\xi = 30\pi" hei 0.6 font texcmr rotate
data res2da
data res2db
d1 lstyle 1 smooth
d2 lstyle 2 smooth
end graph
set font texcmr
set hei 0.7
amove 4.5 28.2
text Sech Pulse with Varying Amplitude (\theta=45^{\circ},\delta=0.5)


There is still (as of 2007) a version of gle available for linux and using the above files I recreated the plot in .eps and convert to a png with

gs -dNOPAUSE -sDEVICE=png256 -sOutputFile=birefring.png -q -dBATCH pulses1.eps

This shows the two birefringent modes traveling forwards and broadening – of course, I didn’t really understand at the time that the oscillations were due to my computational window being too small. But I was already at the limits of what my 386 computer count reasonably handle.