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](http://l.wordpress.com/latex.php?latex=%7B%5Cpartial+u%5Cover%5Cpartial%5Cxi%7D%3D-%5Cdelta%7B%5Cpartial+u%5Cover%5Cpartial%5Ctau%7D%2B%7Bi%5Cover+2%7D%7B%5Cpartial%5E2+u%5Cover%5Cpartial%5Ctau%5E2%7D%2Bi%5Cleft%5B%7Cu%7C%5E2%2B%7B2%5Cover+3%7D%7Cv%7C%5E2%5Cright%5D+u&bg=ffffff&fg=000000&s=3 "{\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](http://l.wordpress.com/latex.php?latex=%7B%5Cpartial+v%5Cover%5Cpartial%5Cxi%7D%3D%2B%5Cdelta%7B%5Cpartial+v%5Cover%5Cpartial%5Ctau%7D%2B%7Bi%5Cover+2%7D%7B%5Cpartial%5E2+v%5Cover%5Cpartial%5Ctau%5E2%7D%2Bi%5Cleft%5B%7Cv%7C%5E2%2B%7B2%5Cover+3%7D%7Cu%7C%5E2%5Cright%5D+v&bg=ffffff&fg=000000&s=3 "{\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.
