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

Birefringent Fibre

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.

Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

%d bloggers like this: