sspropvc.c

用MATLAB仿真光纤的非线性效应中的自相位调制曲线.

/*  File:           sspropvc.c
 *  Authors:        Afrouz Azari ([e-mail])
 *                  Ross A. Pleban ([e-mail])
 *                  Reza Salem ([e-mail])
 *                  Thomas E. Murphy ([e-mail])
 *
 *  Created:        1/17/2001
 *  Modified:       8/18/2006
 *  Version:        3.0
 *  Description:    This file solves the coupled nonlinear
 *                  Schrodinger equations for propagation in an
 *                  optical fiber, using the split-step Fourier
 *                  method.  The routine is compiled as a Matlab
 *                  MEX function that can be invoked directly
 *                  from Matlab.  
 */


/* 
 * USAGE:
 * [u1x,u1y] = sspropvc(u0x,u0y,dt,dz,nz,alphaa,alphab,betapa,betapb,gamma);
 * [u1x,u1y] = sspropvc(u0x,u0y,dt,dz,nz,alphaa,alphab,betapa,betapb,gamma,psp);
 * [u1x,u1y] = sspropvc(u0x,u0y,dt,dz,nz,alphaa,alphab,betapa,betapb,gamma,psp,method);
 * [u1x,u1y] = sspropvc(u0x,u0y,dt,dz,nz,alphaa,alphab,betapa,betapb,gamma,psp,method,maxiter;
 * [u1x,u1y] = sspropvc(u0x,u0y,dt,dz,nz,alphaa,alphab,betapa,betapb,gamma,psp,method,maxiter,tol);
 * sspropvc -option
 *
 * OPTIONS:   (i.e. sspropvc -savewisdom )
 *  -savewisdom
 *  -forgetwisdom
 *  -loadwisdom
 *  -patient
 *  -exhaustive
 *  -measure
 *  -estimate
 */


/*****************************************************************

    Copyright 2006, Thomas E. Murphy

    This file is part of SSPROP.

    SSPROP is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License
    as published by the Free Software Foundation; either version
    2 of the License, or (at your option) any later version.

    SSPROP is distributed in the hope that it will be useful, but
    WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public
    License along with SSPROP; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
    02111-1307 USA

*****************************************************************/

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fftw3.h"
#include "mex.h"

#ifdef SINGLEPREC

#define REAL float
#define COMPLEX fftwf_complex
#define PLAN fftwf_plan
#define MAKE_PLAN fftwf_plan_dft_1d
#define DESTROY_PLAN fftwf_destroy_plan
#define EXECUTE fftwf_execute
#define IMPORT_WISDOM fftwf_import_wisdom_from_file
#define EXPORT_WISDOM fftwf_export_wisdom_to_file
#define FORGET_WISDOM fftwf_forget_wisdom
#define WISFILENAME "fftwf-wisdom.dat"

#else

#define REAL double
#define COMPLEX fftw_complex
#define PLAN fftw_plan
#define MAKE_PLAN fftw_plan_dft_1d
#define DESTROY_PLAN fftw_destroy_plan
#define EXECUTE fftw_execute
#define IMPORT_WISDOM fftw_import_wisdom_from_file
#define EXPORT_WISDOM fftw_export_wisdom_to_file
#define FORGET_WISDOM fftw_forget_wisdom
#define WISFILENAME "fftw-wisdom.dat"

#endif

#define abs2(x) ((*x)[0] * (*x)[0] + (*x)[1] * (*x)[1])
#define round(x) ((int)(x+0.5))
#define pi 3.1415926535897932384626433832795028841972

static int firstcall = 1;       /* =1 when sspropvc first invoked */
int allocated = 0;              /* =1 when memory is allocated */
static int method = FFTW_PATIENT;	/* planner method */

void sspropvc_save_wisdom();
void sspropvc_load_wisdom();
void cscale(COMPLEX*,COMPLEX*,COMPLEX*,COMPLEX*,REAL,int);
void rotate_coord(COMPLEX*,COMPLEX*,const mxArray*,const mxArray*,REAL,REAL,int);
void compute_w(REAL*,REAL,int);
void compute_hahb(COMPLEX*,COMPLEX*,const mxArray*,const mxArray*,const mxArray*,
                  const mxArray*,REAL*,REAL,int);
void compute_H(COMPLEX*,COMPLEX*,COMPLEX*,COMPLEX*,COMPLEX*,COMPLEX*,
               REAL,REAL,int);
void prop_linear_ellipt(COMPLEX* uZa, COMPLEX* uZb, COMPLEX* ha,
                        COMPLEX* hb,COMPLEX* u0a,COMPLEX* u0b,int nt);
void prop_linear_circ(COMPLEX*,COMPLEX*,COMPLEX*,COMPLEX*,COMPLEX*,
                      COMPLEX*,COMPLEX*,COMPLEX*,int);
void nonlinear_propagate(COMPLEX*,COMPLEX*,COMPLEX*,COMPLEX*,COMPLEX*,
                         COMPLEX*,COMPLEX*,COMPLEX*,REAL,REAL,REAL,int);
int is_converged(COMPLEX*,COMPLEX*,COMPLEX*,COMPLEX*,REAL,int);
void inv_rotate_coord(mxArray*,mxArray*,COMPLEX*,COMPLEX*,
                      REAL,REAL,int);
void mexFunction(int, mxArray* [], int, const mxArray* []);


/* Saves the wisdom data to a file specified by WISFILENAME */
void sspropvc_save_wisdom() 
{ 
  FILE *wisfile;
  
  wisfile = fopen(WISFILENAME, "w");
  if (wisfile) {
    mexPrintf("Exporting FFTW wisdom (file = %s).\n", WISFILENAME);
    EXPORT_WISDOM(wisfile);
    fclose(wisfile);
  } 
}


/* Loads the wisdom file into memory specified by WISFILENAME.  The 
 * program automatically attempts to load the wisdom on the first call */
void sspropvc_load_wisdom() 
{ 
  FILE *wisfile;
  
  wisfile = fopen(WISFILENAME, "r");
  if (wisfile) {
	mexPrintf("Importing FFTW wisdom (file = %s).\n", WISFILENAME);
	if (!IMPORT_WISDOM(wisfile)) {
      fclose(wisfile);
      mexErrMsgTxt("could not import wisdom.");
    }
    else
	  fclose(wisfile);
  }
}


/* Assigns a = factor*b  and  x = factor*y  for length nt vectors */
void cscale(COMPLEX* a,COMPLEX* b,COMPLEX* x,COMPLEX* y,REAL factor,int nt)
{
  int jj;
  
  for (jj = 0; jj < nt; jj++) {
    a[jj][0] = factor*b[jj][0];
    a[jj][1] = factor*b[jj][1];
    x[jj][0] = factor*y[jj][0];
    x[jj][1] = factor*y[jj][1];
  }
}


/* Rotates input to the coordinate system defined by chi & psi 
 *
 * Elliptical MATLAB equivalent:
 *   u0a = ( cos(psi)*cos(chi) - j*sin(psi)*sin(chi))*u0x + ...
 *         ( sin(psi)*cos(chi) + j*cos(psi)*sin(chi))*u0y;
 *   u0b = (-sin(psi)*cos(chi) + j*cos(psi)*sin(chi))*u0x + ...
 *         ( cos(psi)*cos(chi) + j*sin(psi)*sin(chi))*u0y;
 * 
 * Circular MATLAB equivalent (when chi = pi/4 and psi = 0):
 *   u0a = (1/sqrt(2)).*(u0x + j*u0y);
 *   u0b = (1/sqrt(2)).*(j*u0x + u0y);
 */
void rotate_coord(COMPLEX* u0a, COMPLEX* u0b,const mxArray* ux, const mxArray* uy,
                  REAL chi, REAL psi, int nt) 
{ 
  REAL cc = (REAL) (cos(psi)*cos(chi));
  REAL ss = (REAL) (sin(psi)*sin(chi));
  REAL sc = (REAL) (sin(psi)*cos(chi));
  REAL cs = (REAL) (cos(psi)*sin(chi));
  int jj;
  if (mxIsComplex(ux) && mxIsComplex(uy))
    for(jj = 0; jj < nt; jj++) {
      u0a[jj][0] = cc*mxGetPr(ux)[jj] + ss* mxGetPi(ux)[jj] + 
                   sc*mxGetPr(uy)[jj] - cs*mxGetPi(uy)[jj];
      u0a[jj][1] = cc*mxGetPi(ux)[jj] - ss*mxGetPr(ux)[jj] + 
                   sc*mxGetPi(uy)[jj] + cs*mxGetPr(uy)[jj];
      u0b[jj][0] = -sc*mxGetPr(ux)[jj] - cs*mxGetPi(ux)[jj] + 
                   cc*mxGetPr(uy)[jj] - ss*mxGetPi(uy)[jj];
      u0b[jj][1] = -sc*mxGetPi(ux)[jj] + cs*mxGetPr(ux)[jj] + 
                   cc*mxGetPi(uy)[jj] + ss*mxGetPr(uy)[jj];
    }
  else if (mxIsComplex(ux))
    for(jj = 0; jj < nt; jj++) {
      u0a[jj][0] = cc*mxGetPr(ux)[jj] + ss* mxGetPi(ux)[jj] + 
                   sc*mxGetPr(uy)[jj];
      u0a[jj][1] = cc*mxGetPi(ux)[jj] - ss*mxGetPr(ux)[jj] + 
                   cs*mxGetPr(uy)[jj];
      u0b[jj][0] = -sc*mxGetPr(ux)[jj] - cs*mxGetPi(ux)[jj] + 
                   cc*mxGetPr(uy)[jj];
      u0b[jj][1] = -sc*mxGetPi(ux)[jj] + cs*mxGetPr(ux)[jj] + 
                   ss*mxGetPr(uy)[jj];
    }
  else if (mxIsComplex(uy))
    for(jj = 0; jj < nt; jj++) {
      u0a[jj][0] = cc*mxGetPr(ux)[jj] +  
                   sc*mxGetPr(uy)[jj] - cs*mxGetPi(uy)[jj];
      u0a[jj][1] = - ss*mxGetPr(ux)[jj] + 
                   sc*mxGetPi(uy)[jj] + cs*mxGetPr(uy)[jj];
      u0b[jj][0] = -sc*mxGetPr(ux)[jj] + 
                   cc*mxGetPr(uy)[jj] - ss*mxGetPi(uy)[jj];
      u0b[jj][1] = cs*mxGetPr(ux)[jj] + 
                   cc*mxGetPi(uy)[jj] + ss*mxGetPr(uy)[jj];
    }
  else 
    for(jj = 0; jj < nt; jj++) {
      u0a[jj][0] = cc*mxGetPr(ux)[jj] +  
                   sc*mxGetPr(uy)[jj];
      u0a[jj][1] = - ss*mxGetPr(ux)[jj] + 
                   cs*mxGetPr(uy)[jj];
      u0b[jj][0] = -sc*mxGetPr(ux)[jj] + 
                   cc*mxGetPr(uy)[jj];
      u0b[jj][1] = cs*mxGetPr(ux)[jj] + 
                   ss*mxGetPr(uy)[jj];
    }
}


/* Compute vector of angular frequency components
 * MATLAB equivalent:  w = wspace(tv); */
void compute_w(REAL* w,REAL dt,int nt)
{
  int jj;
  for (jj = 0; jj <= (nt-1)/2; jj++) {
    w[jj] = 2*pi*jj/(dt*nt);
  }
  for (; jj < nt; jj++) {
    w[jj] = 2*pi*jj/(dt*nt) - 2*pi/dt;
  }
}


/* Compute ha & hb
 * ha = exp[(-alphaa(w)/2 - j*betaa(w))*dz/2])
 * hb = exp[(-alphab(w)/2 - j*betab(w))*dz/2]) 
 *
 * MATLAB Equivalent of ha (similar for hb):
 *   if (length(alphaa) == nt)   % If the user manually specifies alpha(w)
 *     ha = -alphaa/2;
 *   else
 *     ha = 0;
 *     for ii = 0:length(alphaa)-1;
 *       ha = ha - alphaa(ii+1)*(w).^ii/factorial(ii);
 *     end
 *     ha = ha/2;
 *   end
 *
 *   if (length(betapa) == nt)   % If the user manually specifies beta(w)
 *     ha = ha - j*betapa;
 *   else
 *     for ii = 0:length(betapa)-1;
 *       ha = ha - j*betapa(ii+1)*(w).^ii/factorial(ii);
 *     end
 *   end
 *   ha = exp(ha.*dz/2); % ha = exp[(-alphaa/2 - j*betaa)*dz/2])
 */
void compute_hahb(COMPLEX* ha,COMPLEX* hb,const mxArray* mxAlphaa,
                  const mxArray* mxAlphab,const mxArray* mxBetaa,const mxArray* mxBetab,
                  REAL* w,REAL dz,int nt)
{
  int nalphaa,nalphab,nbetaa,nbetab;    /* # of elements */
  double *alphaa,*alphab,*betaa,*betab; /* taylor coefficients */
  REAL fii,wii,aa,ab,phasea,phaseb;     /* temporary variables */
  int jj,ii;                            /* counters */
  
  nalphaa = mxGetNumberOfElements(mxAlphaa);
  nalphab = mxGetNumberOfElements(mxAlphab);
  nbetaa = mxGetNumberOfElements(mxBetaa);
  nbetab = mxGetNumberOfElements(mxBetab);
  alphaa = mxGetPr(mxAlphaa);
  alphab = mxGetPr(mxAlphab);
  betaa = mxGetPr(mxBetaa);
  betab = mxGetPr(mxBetab);
  
  for (jj = 0; jj < nt; jj++) {
    if (nalphaa != nt)
	  for (ii = 0, aa = 0, fii = 1, wii = 1; 
		   ii < nalphaa; 
		   ii++, fii*=ii, wii*=w[jj]) 
		aa += wii*((REAL) alphaa[ii])/fii;
  	else
	  aa = (REAL)alphaa[jj];
    if (nalphab != nt)
	  for (ii = 0, ab = 0, fii = 1, wii = 1; 
		   ii < nalphab; 
		   ii++, fii*=ii, wii*=w[jj]) 
		ab += wii*((REAL) alphab[ii])/fii;
  	else
	  ab = (REAL)alphab[jj];
    if (nbetaa != nt) 	 
	  for (ii = 0, phasea = 0, fii = 1, wii = 1; 
		   ii < nbetaa; 
		   ii++, fii*=ii, wii*=w[jj]) 
		phasea += wii*((REAL)betaa[ii])/fii;
  	else 
	  phasea = (REAL)betaa[jj];
    if (nbetab != nt) 	 
	  for (ii = 0, phaseb = 0, fii = 1, wii = 1; 
		   ii < nbetab; 
		   ii++, fii*=ii, wii*=w[jj]) 
		phaseb += wii*((REAL)betab[ii])/fii;
  	else 
	  phaseb = (REAL)betab[jj];
    ha[jj][0] = +exp(-aa*dz/4)*cos(phasea*dz/2);
    ha[jj][1] = -exp(-aa*dz/4)*sin(phasea*dz/2);
    hb[jj][0] = +exp(-ab*dz/4)*cos(phaseb*dz/2);
    hb[jj][1] = -exp(-ab*dz/4)*sin(phaseb*dz/2);
  }
}


/* Compute H matrix = [ h11 h12 
 *                      h21 h22 ] for linear propagation 
 *
 * MATLAB Equivalent:
 *   h11 = ( (1+sin(2*chi))*ha + (1-sin(2*chi))*hb )/2;
 *   h12 = -j*exp(+j*2*psi)*cos(2*chi)*(ha-hb)/2;
 *   h21 = +j*exp(-j*2*psi)*cos(2*chi)*(ha-hb)/2;
 *   h22 = ( (1-sin(2*chi))*ha + (1+sin(2*chi))*hb )/2;
 */
void compute_H(COMPLEX* h11,COMPLEX* h12,COMPLEX* h21,COMPLEX* h22,
               COMPLEX* ha,COMPLEX* hb,REAL chi,REAL psi,int nt)
{
  int jj;
  REAL halfPsin,halfMsin,sincos,coscos;
  halfPsin = .5 + .5*sin(2*chi);
  halfMsin = .5 - .5*sin(2*chi);
  sincos = .5*sin(2*psi)*cos(2*chi);
  coscos = .5*cos(2*psi)*cos(2*chi);
  
  for(jj = 0; jj < nt; jj++)
  {
    h11[jj][0] = halfPsin*ha[jj][0] + halfMsin*hb[jj][0];
    h11[jj][1] = halfPsin*ha[jj][1] + halfMsin*hb[jj][1];
    h12[jj][0] = sincos*(ha[jj][0]-hb[jj][0])+ 
                 coscos*(ha[jj][1]-hb[jj][1]);
    h12[jj][1] = sincos*(ha[jj][1]-hb[jj][1])- 
                 coscos*(ha[jj][0]-hb[jj][0]);
    h21[jj][0] = sincos*(ha[jj][0]-hb[jj][0]) -
                 coscos*(ha[jj][1]-hb[jj][1]);
    h21[jj][1] = sincos*(ha[jj][1]-hb[jj][1]) +
                 coscos*(ha[jj][0]-hb[jj][0]);
    h22[jj][0] = halfMsin*ha[jj][0] + halfPsin*hb[jj][0];
    h22[jj][1] = halfMsin*ha[jj][1] + halfPsin*hb[jj][1];
  }
}


/* Computes elliptical linear progation according to 
 * the matrix multiplication of:
 *   [ Ua(Z) ] = [ ha  0  ] * [ Ua(0) ]
 *   [ Ub(Z) ]   [ 0   hb ]   [ Ub(0) ] 
 *
 * MATLAB Equivalent:
 *   uZa = ha .* u0a;
 *   uZb = hb .* u0b;
 */
void prop_linear_ellipt(COMPLEX* uZa, COMPLEX* uZb, COMPLEX* ha,
                        COMPLEX* hb,COMPLEX* u0a,COMPLEX* u0b,int nt)
{
  int jj;

  for (jj = 0; jj < nt; jj++) {
    uZa[jj][0] = ha[jj][0]*u0a[jj][0] - ha[jj][1]*u0a[jj][1] ;
    uZa[jj][1] = ha[jj][0]*u0a[jj][1] + ha[jj][1]*u0a[jj][0] ;
    uZb[jj][0] = hb[jj][0]*u0b[jj][0] - hb[jj][1]*u0b[jj][1] ;
    uZb[jj][1] = hb[jj][0]*u0b[jj][1] + hb[jj][1]*u0b[jj][0] ;
  }
}


/* Computes the circular linear progation according to the 
 * matrix multiplication of:
 *   [ Ua(Z) ] = [ h11 h12 ] * [ Ua(0) ]
 *   [ Ub(Z) ] = [ h21 h22 ]   [ Ub(0) ] 
 *
 * MATLAB Equivalent:
 *   uZa = h11 .* u0a + h12 .* u0b;
 *   uZb = h21 .* u0a + h22 .* u0b;
 */
void prop_linear_circ(COMPLEX* uZa, COMPLEX* uZb, COMPLEX* h11,
                      COMPLEX* h12,COMPLEX* h21,COMPLEX* h22,COMPLEX* u0a,
                      COMPLEX* u0b,int nt)
{
  int jj;

  for (jj = 0; jj < nt; jj++) {
    uZa[jj][0] = h11[jj][0]*u0a[jj][0] + h12[jj][0]*u0b[jj][0] -
                 h11[jj][1]*u0a[jj][1] - h12[jj][1]*u0b[jj][1];
    uZa[jj][1] = h11[jj][0]*u0a[jj][1] + h12[jj][0]*u0b[jj][1] +
                 h11[jj][1]*u0a[jj][0] + h12[jj][1]*u0b[jj][0];
    uZb[jj][0] = h21[jj][0]*u0a[jj][0] + h22[jj][0]*u0b[jj][0] -
                 h21[jj][1]*u0a[jj][1] - h22[jj][1]*u0b[jj][1];
    uZb[jj][1] = h21[jj][0]*u0a[jj][1] + h22[jj][0]*u0b[jj][1] +
                 h21[jj][1]*u0a[jj][0] + h22[jj][1]*u0b[jj][0];
  }
}


/* Computes nonlinear propagation according to the following equations:
 *
 * Elliptical Equivalent:
 * dua/dz = (-j*gamma/3)*[(2+cos(2X)^2*|ua|^2 + (2+2sin(2X)^2)*|ub|^2] * ua
 * dub/dz = (-j*gamma/3)*[(2+cos(2X)^2*|ub|^2 + (2+2sin(2X)^2)*|ua|^2] * ub
 *
 * Circular Equivalent:
 * dua/dz = (-j*2*gamma/3)*(|ua|^2 + 2*|ub|^2)*ua
 * dub/dz = (-j*2*gamma/3)*(|ub|^2 + 2*|ua|^2)*ub
 */
void nonlinear_propagate(COMPLEX* uva,COMPLEX* uvb,COMPLEX* uahalf,
                    COMPLEX* ubhalf,COMPLEX* u0a,COMPLEX* u0b,
                    COMPLEX* u1a,COMPLEX* u1b,REAL gamma,REAL dz,
                    REAL chi, int nt)
{
  int jj;
  REAL coef,twoPcos,twoPsin;
  coef = (REAL) ((1.0/3.0)*gamma*dz);