sspropvc.c

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

  twoPcos = (2 + cos(2*chi)*cos(2*chi)) / 2;
  twoPsin = (2 + 2*sin(2*chi)*sin(2*chi)) / 2;
    
  for(jj = 0; jj < nt; jj++) {
    uva[jj][0] = uahalf[jj][0]*cos(coef*(
                   twoPcos*(abs2(&u0a[jj])+abs2(&u1a[jj])) +
                   twoPsin*(abs2(&u0b[jj])+abs2(&u1b[jj])) )) 
               + uahalf[jj][1]*sin(coef*(
                   twoPcos*(abs2(&u0a[jj])+abs2(&u1a[jj])) +
                   twoPsin*(abs2(&u0b[jj])+abs2(&u1b[jj])) ));
    uva[jj][1] = uahalf[jj][1]*cos(coef*(
                   twoPcos*(abs2(&u0a[jj])+abs2(&u1a[jj])) +
                   twoPsin*(abs2(&u0b[jj])+abs2(&u1b[jj])) )) 
               - uahalf[jj][0]*sin(coef*(
                   twoPcos*(abs2(&u0a[jj])+abs2(&u1a[jj])) +
                   twoPsin*(abs2(&u0b[jj])+abs2(&u1b[jj])) ));
    uvb[jj][0] = ubhalf[jj][0]*cos(coef*(
                   twoPcos*(abs2(&u0b[jj])+abs2(&u1b[jj])) +
                   twoPsin*(abs2(&u0a[jj])+abs2(&u1a[jj])) )) 
               + ubhalf[jj][1]*sin(coef*(
                   twoPcos*(abs2(&u0b[jj])+abs2(&u1b[jj])) +
                   twoPsin*(abs2(&u0a[jj])+abs2(&u1a[jj])) ));
    uvb[jj][1] = ubhalf[jj][1]*cos(coef*(
                   twoPcos*(abs2(&u0b[jj])+abs2(&u1b[jj])) +
                   twoPsin*(abs2(&u0a[jj])+abs2(&u1a[jj])) )) 
               - ubhalf[jj][0]*sin(coef*(
                   twoPcos*(abs2(&u0b[jj])+abs2(&u1b[jj])) +
                   twoPsin*(abs2(&u0a[jj])+abs2(&u1a[jj])) ));
  }
}


/* Returns non-zero if uva & uvb have converged towards u1a & u1b with
 * a tolerance less than tol 
 *
 * MATLAB equivalent:
 *   ( sqrt(norm(uva-u1a,2).^2+norm(uvb-u1b,2).^2) / ...
 *      sqrt(norm(u1a,2).^2+norm(u1b,2).^2) ) < tol 
 */
int is_converged(COMPLEX* uva,COMPLEX* u1a,COMPLEX* uvb,COMPLEX* u1b,
                 REAL tol,int nt) 
{ 
  int jj;
  REAL num,denom;
  for(jj = 0, num = 0, denom = 0; jj < nt; jj++) {
    num += (uva[jj][0]/nt-u1a[jj][0])*(uva[jj][0]/nt-u1a[jj][0]) +  
           (uva[jj][1]/nt-u1a[jj][1])*(uva[jj][1]/nt-u1a[jj][1]) +
           (uvb[jj][0]/nt-u1b[jj][0])*(uvb[jj][0]/nt-u1b[jj][0]) +  
           (uvb[jj][1]/nt-u1b[jj][1])*(uvb[jj][1]/nt-u1b[jj][1]);
    denom += abs2(&u1a[jj]) + abs2(&u1b[jj]);
  }
  return ( sqrt(num)/sqrt(denom) < tol);
}


/* Rotates back to input coordinate system, where u1x & u1y are the 
 * outputs and u1a and u1b are the inputs 
 *
 * Elliptical equivalent:
 *   u1x = ( cos(psi)*cos(chi) + j*sin(psi)*sin(chi))*u1a + ...
 *         (-sin(psi)*cos(chi) - j*cos(psi)*sin(chi))*u1b;
 *   u1y = ( sin(psi)*cos(chi) - j*cos(psi)*sin(chi))*u1a + ...
 *         ( cos(psi)*cos(chi) - j*sin(psi)*sin(chi))*u1b;
 * 
 * Circular MATLAB equivalent (when chi = pi/4 and psi = 0):
 *   u1x = (1/sqrt(2)).*(u1a-j*u1b) ;
 *   u1y = (1/sqrt(2)).*(-j*u1a+u1b) ;
 */
void inv_rotate_coord(mxArray* u1x, mxArray* u1y,COMPLEX* u1a,COMPLEX* u1b,
                      REAL chi, REAL psi, int nt) 
{ 
  REAL cc = cos(psi)*cos(chi);
  REAL ss = sin(psi)*sin(chi);
  REAL sc = sin(psi)*cos(chi);
  REAL cs = cos(psi)*sin(chi);
  int jj;
  for(jj = 0; jj < nt; jj++) {
    mxGetPr(u1x)[jj] = cc*u1a[jj][0] - ss*u1a[jj][1] -
                       sc*u1b[jj][0] + cs*u1b[jj][1];
    mxGetPi(u1x)[jj] = cc*u1a[jj][1] + ss*u1a[jj][0] -
                       sc*u1b[jj][1] - cs*u1b[jj][0];
    mxGetPr(u1y)[jj] = sc*u1a[jj][0] + cs*u1a[jj][1] +
                       cc*u1b[jj][0] + ss*u1b[jj][1];
    mxGetPi(u1y)[jj] = sc*u1a[jj][1] - cs*u1a[jj][0] +
                       cc*u1b[jj][1] - ss*u1b[jj][0];
  }
}


/* This is the gateway function between MATLAB and SSPROPVC.  It
 * serves as the main(). */
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{ 
  COMPLEX *u0a, *u0b, *uafft, *ubfft, *uahalf, *ubhalf,
          *uva, *uvb, *u1a, *u1b;
  
  COMPLEX *ha, *hb;  /* exp{ (-Alpha(w)/2-jBeta(w)) z} */
  COMPLEX *h11, *h12,/* linear propgation coefficients */
          *h21, *h22;
    
  REAL dt;           /* time step */
  REAL dz;           /* propagation stepsize */
  int nz;            /* number of z steps to take */
  REAL gamma;        /* nonlinearity coefficient */
  REAL chi = 0.0;    /* degree of ellipticity  */
  REAL psi = 0.0;    /* angular orientation to x-axis  */
  int maxiter = 4;   /* max number of iterations */
  REAL tol = 1e-5;   /* convergence tolerance */

  int nt;            /* number of fft points */
  
  REAL* w;           /* vector of angular frequencies */

  PLAN p1a,p1b,ip1a,ip1b;   /* fft plans for 1st linear half */
  PLAN p2a,p2b,ip2a,ip2b;   /* fft plans for 2nd linear half */
  
  int converged;            /* holds the return of is_converged */
  char methodstr[11];       /* method name: 'circular or 'elliptical' */
  int elliptical = 1;       /* if elliptical method, then != 0 */

  char argstr[100];	 /* string argument */
  
  int iz,ii,jj;      /* loop counters */
  
  if (nrhs == 1) {
	if (mxGetString(prhs[0],argstr,100)) 
	  mexErrMsgTxt("Unrecognized option.");
	
	if (!strcmp(argstr,"-savewisdom")) {
	  sspropvc_save_wisdom();
	}
	else if (!strcmp(argstr,"-forgetwisdom")) {
	  FORGET_WISDOM();
	}
	else if (!strcmp(argstr,"-loadwisdom")) {
	  sspropvc_load_wisdom();
	}
	else if (!strcmp(argstr,"-patient")) {
	  method = FFTW_PATIENT;
	}
	else if (!strcmp(argstr,"-exhaustive")) {
	  method = FFTW_EXHAUSTIVE;
	}
	else if (!strcmp(argstr,"-measure")) {
	  method = FFTW_MEASURE;
	}
	else if (!strcmp(argstr,"-estimate")) {
	  method = FFTW_ESTIMATE;
	}
	else
	  mexErrMsgTxt("Unrecognized option.");
	return;
  }
  
  if (nrhs < 10) 
    mexErrMsgTxt("Not enough input arguments provided.");
  if (nlhs > 2)
    mexErrMsgTxt("Too many output arguments.");
  
  if (firstcall) {  /* attempt to load wisdom file on first call */
	sspropvc_load_wisdom();
    firstcall = 0;
  }

  /* parse input arguments */
  dt = (REAL) mxGetScalar(prhs[2]);
  dz = (REAL) mxGetScalar(prhs[3]);
  nz = round(mxGetScalar(prhs[4]));
  gamma = (REAL) mxGetScalar(prhs[9]);

  if (nrhs > 10) { /* default is chi = psi = 0.0 */
    psi = (REAL) mxGetScalar(prhs[10]); 
	if (mxGetNumberOfElements(prhs[10]) > 1)
	  chi = (REAL) (mxGetPr(prhs[10])[1]); 
  } 
 
  if (nrhs > 11) { /* default method is elliptical */
    if (mxGetString(prhs[11],methodstr,11)) /* fail */
      mexErrMsgTxt("incorrect method: elliptical or ciruclar only");
    else { /* success */
      if (!strcmp(methodstr,"circular"))
        elliptical = 0;
      else if(!strcmp(methodstr,"elliptical"))
        elliptical = 1;
      else
         mexErrMsgTxt("incorrect method: elliptical or ciruclar only");
    }
  }
    
  if (nrhs > 12) /* default = 4 */
	maxiter = round(mxGetScalar(prhs[12]));
  
  if (nrhs > 13) /* default = 1e-5 */
	tol = (REAL) mxGetScalar(prhs[13]);

  nt = mxGetNumberOfElements(prhs[0]);  /* # of points in vectors */
  
  /* allocate memory */
  u0a = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  u0b = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  uafft = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  ubfft = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  uahalf = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  ubhalf = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  uva = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  uvb = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  u1a = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  u1b = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  ha = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  hb = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  h11 = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  h12 = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  h21 = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  h22 = (COMPLEX*) mxMalloc(sizeof(COMPLEX)*nt);
  w = (REAL*)mxMalloc(sizeof(REAL)*nt);
  plhs[0] = mxCreateDoubleMatrix(nt,1,mxCOMPLEX);
  plhs[1] = mxCreateDoubleMatrix(nt,1,mxCOMPLEX);
  
  /* fftw3 plans */
  p1a = MAKE_PLAN(nt, u0a, uafft, FFTW_FORWARD, method);
  p1b = MAKE_PLAN(nt, u0b, ubfft, FFTW_FORWARD, method);
  ip1a = MAKE_PLAN(nt, uahalf, uahalf, FFTW_BACKWARD, method);
  ip1b = MAKE_PLAN(nt, ubhalf, ubhalf, FFTW_BACKWARD, method);
  p2a = MAKE_PLAN(nt, uva, uva, FFTW_FORWARD, method);
  p2b = MAKE_PLAN(nt, uvb, uvb, FFTW_FORWARD, method);
  ip2a = MAKE_PLAN(nt, uafft, uva, FFTW_BACKWARD, method);
  ip2b = MAKE_PLAN(nt, ubfft, uvb, FFTW_BACKWARD, method);

  allocated = 1;
  
  /* Compute vector of angular frequency components
   * MATLAB equivalent:  w = wspace(tv); */
  compute_w(w,dt,nt);
  
  /* Compute ha & hb vectors
   * ha = exp[(-alphaa(w)/2 - j*betaa(w))*dz/2])
   * hb = exp[(-alphab(w)/2 - j*betab(w))*dz/2]) 
   * prhs[5]=alphaa  prhs[6]=alphab  prhs[7]=betaa  prhs[8]=betab */
  compute_hahb(ha,hb,prhs[5],prhs[6],prhs[7],prhs[8],w,dz,nt);
  
  mexPrintf("Performing split-step iterations ... ");
  
  if (elliptical) { /* Elliptical Method */
    
    /* Rotate to eignestates of fiber 
     *   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;
     */
    rotate_coord(u0a,u0b,prhs[0],prhs[1],chi,psi,nt);
      
    cscale(u1a,u0a,u1b,u0b,1.0,nt); /* u1a=u0a  u1b=u0b */
    
    EXECUTE(p1a);  /* uafft = fft(u0a) */
    EXECUTE(p1b);  /* ubfft = fft(u0b) */
    
    for(iz=1; iz <= nz; iz++)
    {
      /* Linear propagation (1st half):
       * uahalf = ha .* uafft
       * ubhalf = hb .* ubfft */
      prop_linear_ellipt(uahalf,ubhalf,ha,hb,uafft,ubfft,nt);
      
      EXECUTE(ip1a);  /* uahalf = ifft(uahalf) */
      EXECUTE(ip1b);  /* ubhalf = ifft(ubhalf) */
      
      /* uahalf=uahalf/nt  ubhalf=ubhalf/nt */
      cscale(uahalf,uahalf,ubhalf,ubhalf,1.0/nt,nt);
  
      ii = 0;
      do
      {
        /* Calculate nonlinear section: output=uva,uvb */
        nonlinear_propagate(uva,uvb,uahalf,ubhalf,u0a,u0b,u1a,u1b,
                            gamma,dz,chi,nt);
        
      
        EXECUTE(p2a);  /* uva = fft(uva) */
        EXECUTE(p2b);  /* uvb = fft(uvb) */
      
        /* Linear propagation (2nd half):
         * uafft = ha .* uva
         * ubfft = hb .* uvb */
         prop_linear_ellipt(uafft,ubfft,ha,hb,uva,uvb,nt);
     
        EXECUTE(ip2a);  /* uva = ifft(uafft) */
        EXECUTE(ip2b);  /* uvb = ifft(ubfft) */
        
        /* Check if uva & u1a  and  uvb & u1b converged 
         * converged = ( ( sqrt(norm(uva-u1a,2).^2+norm(uvb-u1b,2).^2) /...
         *                 sqrt(norm(u1a,2).^2+norm(u1b,2).^2) ) < tol )
         */
        converged = is_converged(uva,u1a,uvb,u1b,tol,nt);
      
        /* u1a=uva/nt  u1b=uvb/nt */
        cscale(u1a,uva,u1b,uvb,1.0/nt,nt);
      
        ii++;
      } while(!converged && ii < maxiter);  /* end convergence loop */
    
      if(ii == maxiter)
        mexPrintf("Warning: Failed to converge to %f in %d iterations\n",
                  tol,maxiter);
    
      /* u0a=u1a  u0b=u1b */
      cscale(u0a,u1a,u0b,u1b,1.0,nt);

    } /* end step loop */
    
    /* Rotate back to original x-y basis
     *  u1x = ( cos(psi)*cos(chi) + j*sin(psi)*sin(chi))*u1a + ...
     *        (-sin(psi)*cos(chi) - j*cos(psi)*sin(chi))*u1b;
     *  u1y = ( sin(psi)*cos(chi) - j*cos(psi)*sin(chi))*u1a + ...
     *        ( cos(psi)*cos(chi) - j*sin(psi)*sin(chi))*u1b;
     */
    inv_rotate_coord(plhs[0],plhs[1],u1a,u1b,chi,psi,nt);
    
  } 
  else {  /* Circular method */ 
    
    /* Compute H matrix = [ h11 h12 
     *                      h21 h22 ] for linear propagation
     *   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;
     */
    compute_H(h11,h12,h21,h22,ha,hb,chi,psi,nt);
      
    /* Rotate to circular coordinate system 
     *   u0a = (1/sqrt(2)).*(u0x + j*u0y);
     *   u0b = (1/sqrt(2)).*(j*u0x + u0y); */
    rotate_coord(u0a,u0b,prhs[0],prhs[1],pi/4,0,nt);
      
    cscale(u1a,u0a,u1b,u0b,1.0,nt); /* u1a=u0a  u1b=u0b */
    
    EXECUTE(p1a);  /* uafft = fft(u0a) */
    EXECUTE(p1b);  /* ubfft = fft(u0b) */
      
    for(iz=1; iz <= nz; iz++)
    {
      /* Linear propagation (1st half):
       * uahalf = h11 .* uafft + h12 .* ubfft
       * ubhalf = h21 .* uafft + h22 .* ubfft */
      prop_linear_circ(uahalf,ubhalf,h11,h12,h21,h22,uafft,ubfft,nt);
      
      EXECUTE(ip1a);  /* uahalf = ifft(uahalf) */
      EXECUTE(ip1b);  /* ubhalf = ifft(ubhalf) */
      
      /* uahalf=uahalf/nt  ubhalf=ubhalf/nt */
      cscale(uahalf,uahalf,ubhalf,ubhalf,1.0/nt,nt);
  
      ii = 0;
      do
      {
        /* Calculate nonlinear section: output=uva,uvb */
         nonlinear_propagate(uva,uvb,uahalf,ubhalf,u0a,u0b,u1a,u1b,
                             gamma,dz,pi/4,nt);
      
        EXECUTE(p2a);  /* uva = fft(uva) */
        EXECUTE(p2b);  /* uvb = fft(uvb) */
      
        /* Linear propagation (2nd half):
         * uafft = h11 .* uva + h12 .* uvb
         * ubfft = h21 .* uva + h22 .* uvb */
        prop_linear_circ(uafft,ubfft,h11,h12,h21,h22,uva,uvb,nt);
     
        EXECUTE(ip2a);  /* uva = ifft(uafft) */
        EXECUTE(ip2b);  /* uvb = ifft(ubfft) */
      
        /* Check if uva & u1a  and  uvb & u1b converged 
         *   ( sqrt(norm(uva-u1a,2).^2+norm(uvb-u1b,2).^2) /...
         *     sqrt(norm(u1a,2).^2+norm(u1b,2).^2) ) < tol
         */
        converged = is_converged(uva,u1a,uvb,u1b,tol,nt);
      
        /* u1a=uva/nt  u1b=uvb/nt */
        cscale(u1a,uva,u1b,uvb,1.0/nt,nt);
      
        ii++;
      } while(!converged && ii < maxiter);  /* end convergence loop */
    
      if(ii == maxiter)
        mexPrintf("Warning: Failed to converge to %f in %d iterations\n",
                  tol,maxiter);
    
      /* u0a=u1a  u0b=u1b */
      cscale(u0a,u1a,u0b,u1b,1.0,nt);

    } /* end step loop */
    
    /* Rotate back to orignal x-y basis
     *   u1x = (1/sqrt(2)).*(u1a-j*u1b) ;
     *   u1y = (1/sqrt(2)).*(-j*u1a+u1b) ; */
    inv_rotate_coord(plhs[0],plhs[1],u1a,u1b,pi/4,0,nt);
    
  } /* end circular method */      
  

  mexPrintf("done.\n");

  if (allocated) {
    /* destroy fftw3 plans */
    DESTROY_PLAN(p1a);
    DESTROY_PLAN(p1b);
    DESTROY_PLAN(ip1a);
    DESTROY_PLAN(ip1b);
    DESTROY_PLAN(p2a);
    DESTROY_PLAN(p2b);
    DESTROY_PLAN(ip2a);
    DESTROY_PLAN(ip2b);

    /* de-allocate memory */
    mxFree(u0a);
    mxFree(u0b);
    mxFree(uafft);
    mxFree(ubfft);
    mxFree(uahalf);
    mxFree(ubhalf);
    mxFree(uva);
    mxFree(uvb);
    mxFree(u1a);
    mxFree(u1b);
    mxFree(ha);
    mxFree(hb);
    mxFree(h11);
    mxFree(h12);
    mxFree(h21);
    mxFree(h22);
    mxFree(w);
    
    allocated = 0;
  }
} /* end mexFunction */