Rcjp's Weblog

April 5, 1993

Soliton Program in C

Filed under: c — rcjp @ 7:36 pm

/*
   PROPAGATES TO GENERATE A 3RD ORDER SOLITON USING BEAM PROPAGATION
   METHOD (SPLIT STEP FOURIER METHOD) OVER ONE SOLITON PERIOD
   ALSO SHOWS THE TIME FOR EXECUTION IN TIME.LOG
   COMPILE WITH FULL OPTIMIZATION
   */
#include <time.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>

double u0[600], uz[600], uzh[600], u02[600];

int main()
{
  time_t start, finish;
  double rt, sb2, vwind, pi, pioverrt, rtover2, w, theta;
  double h, coeff2, mytime, c, s, temp;
  int numz, nowrit, nn, n, i, iycount, isign, idummy, j;
  FILE *solitondata, *logfile;

  void fourier(double *data, int nn, int isign);

  if ( (solitondata = fopen("sol3", "wt")) == NULL ){
    fprintf(stderr, "Cannot open data file");
    exit(1);
  }
  if ( (logfile = fopen("time.log", "wt")) == NULL ){
    fprintf(stderr, "Cannot open log file");
    exit(1);
  }
  start = time(NULL);
  numz  = 1;             /* number of slices */
  nn            = 256;   /* number of points across pulse */
  sb2   = -1.0;  /* sign of beta 2 */
  n             = 3;     /* order of soliton */
  rt            = 8.0;   /* width of pulse */
  vwind = 4.0;   /* view window (pulse width of ouputted results) */
  nowrit= 100;   /* nowrit of iterations where no result written */

  pi       = 4.0 * atan(1.0);
  pioverrt = pi / rt;
  rtover2  = rt / 2.0;
  h                     = pi/(2.0 * numz * nowrit);
  coeff2        = sb2 * h / 4.0;

  /* Generate Pulse */
  for(i=1; i<2*nn; i+=2){
    mytime  = rtover2 * (i-nn)/(nn-1);
    u0[i]   = 1.0/cosh(mytime);
    u02[i]  = u0[i]*u0[i];
  }

  /* Copy initial pulse to uz() which is used to calculate the
     pulse shape at z+h/2 and uzh */
  for(i=1; i<=2*nn; i++)
    uzh[i] = uz[i] = u0[i];

  /* Start of main loop - propagate pulse */
  for(iycount=1; iycount<=numz; iycount++){

    /* Write part of pulse inside view window to file */
    for(i=1; i<2*nn; i+=2){
      mytime = rtover2 * (i-nn)/(nn-1);
      if (fabs(mytime) <= vwind/2.0)
        fprintf(solitondata, "%20e%20e\n", mytime, u02[i]);
    }
    printf("Soliton Order %d, %d of %d\n", n, iycount, numz);

    /* Dummy loop propagates pulse without writing out results */
    for(idummy=1; idummy<=nowrit; idummy++){

      /* Calculate dispersive effects between z->z+h/2 */
      isign = 1;
      fourier(uz, nn, isign);
      for(i=1; i<=nn+1; i+=2){
        w           = pioverrt*(i-1);
        theta       = coeff2*w*w;
        c           = cos(theta);
        s           = sin(theta);
        temp        = uz[i];
        uz[i]       = uz[i]*c-uz[i+1]*s;
        uz[i+1]     = temp*s+uz[i+1]*c;
      }
      for(i=nn+3; i<=2*nn-1; i+=2){
        w           = pioverrt*(i-2*nn-1);
        theta       = coeff2*w*w;
        c           = cos(theta);
        s           = sin(theta);
        temp        = uz[i];
        uz[i]       = uz[i]*c-uz[i+1]*s;
        uz[i+1]     = temp*s+uz[i+1]*c;
      }
      isign = -1;
      fourier(uz, nn, isign);

      /* Multiply by 1/nn */
      for(i=1; i<=2*nn; i++)
        uz[i] /= nn;

      /* Calculate nonlinear effects.
         J-Loop used to find the pulse shape at z+h
         starting with the approximation that it is the pulse in z*/
      for(j=1; j<=3; j++){
        for (i=1; i<2*nn; i+=2){
          theta         = 0.5*h*n*n*(u02[i]+uzh[i]*uzh[i]+uzh[i+1]*uzh[i+1]);
          c                     = cos(theta);
          s                     = sin(theta);
          /* Multiply with dispersion only term uz[] */
          uzh[i]        = uz[i]*c-uz[i+1]*s;
          uzh[i+1]      = uz[i]*s+uz[i+1]*c;
        }
        /* Calculate dispersive effects between z+h/2 -> z+h */
        isign = 1;
        fourier(uzh, nn, isign);

        for(i=1; i<=nn+1; i+=2){
          w           = pioverrt*(i-1);
          theta       = coeff2*w*w;
          c           = cos(theta);
          s           = sin(theta);
          temp        = uzh[i];
          uzh[i]      = uzh[i]*c-uzh[i+1]*s;
          uzh[i+1]    = temp*s+uzh[i+1]*c;
        }
        for(i=nn+3; i<=2*nn-1; i+=2){
          w           = pioverrt*(i-2*nn-1);
          theta       = coeff2*w*w;
          c           = cos(theta);
          s           = sin(theta);
          temp        = uzh[i];
          uzh[i]      = uzh[i]*c-uzh[i+1]*s;
          uzh[i+1]    = temp*s+uzh[i+1]*c;
        }
        isign = -1;
        fourier(uzh, nn, isign);

        /* Multiply by 1/nn */
        for(i=1; i<=2*nn; i++)
          uzh[i] /= nn;
      }
      /* Re-initialise u0[], uz[] and uzh[] to be the new pulse shape*/
      for(i=1; i<=2*nn; i++)
        u0[i] = uz[i] = uzh[i];
      for(i=1; i<2*nn; i+=2)
        u02[i] = u0[i]*u0[i]+u0[i+1]*u0[i+1];
    }
  }
  finish = time(NULL);
  fprintf(logfile,"Time for execution of soliton.c was %f seconds\n",
          difftime(finish, start));
  fprintf(logfile,"Number of slices %d\n", numz);
  fprintf(logfile,"Number of dummy loops %d\n", nowrit);
  fprintf(logfile,"Number of points across pulse %d\n", nn);
  fclose(solitondata);
  fclose(logfile);
  return 0;
}


void fourier(double *data, int nn, int isign)
{
  int i, j, istep, m, mmax, n;
  double tempi, tempr, theta, wi, wpi, wpr, wr, wtemp;

  n = nn << 1;
  j = 1;
  for(i=1; i < n; i+=2){
    if(j > i){
      tempr       = data[j];
      tempi       = data[j+1];
      data[j]     = data[i];
      data[j+1]   = data[i+1];
      data[i]     = tempr;
      data[i+1]   = tempi;
    }
    m = n >> 1;
    while(m >= 2 && j > m){
      j -= m;
      m >>= 1;
    }
    j += m;
  }
  mmax = 2;
  while(n > mmax){
    istep   = mmax << 1;
    theta   = 6.28318530717959/(isign*mmax);
    wpr     =-2.0 * pow(sin(0.5*theta),2.0);
    wpi     = sin(theta);
    wr      = 1.0;
    wi      = 0.0;
    for(m = 1; m < mmax; m+=2){
      for(i = m; i <= n; i+=istep){
        j           = i+mmax;
        tempr       = wr*data[j] - wi*data[j+1];
        tempi       = wr*data[j+1] + wi*data[j];
        data[j]     = data[i]-tempr;
        data[j+1]   = data[i+1]-tempi;
        data[i]     += tempr;
        data[i+1]   += tempi;
      }
      wtemp   = wr;
      wr      += wr*wpr-wi*wpi;
      wi      += wi*wpr+wtemp*wpi;
    }
    mmax = istep;
  }
  return;
}

Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post.

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: