Fresnel積分の基本形

Fresnel（フレネル）積分とは
$$
\int_0^{x} e^{i\frac{\pi}{2}\sigma^{2}}\mathrm{d}\sigma=c(x)+is(x)
$$
または
$$
\begin{align}
\int_0^{x}\cos\frac{\pi}{2}\sigma^{2}\mathrm{d}\sigma&=c(x) \
\int_0^{x}\sin\frac{\pi}{2}\sigma^{2}\mathrm{d}\sigma&=s(x)
\end{align}
$$
と書いた時の$c(x)$と$s(x)$の組のことです。これが何を意味するかについてはWikipediaに委ねます。いわゆる特殊関数に含まれているもので、解析的に計算することができません。ただし数値計算法はよく知られており、Numerical Recipesにも書かれています。

参考までに、ZMの中での実装コードを示します。

# define Z_FRESNEL_FP_MIN ( 1.0e-30 )

/* Fresnel integral by Gilbert formula */
static bool _zFresnelIntgGilbert(double x, double *s, double *c)
{
  int i, n;
  double sum, t, t0;

  *s = sum = 0;
  *c = t = x;
  t0 = zPI_2 * x * x;
  for( n=3, i=1; i<=Z_MAX_ITER_NUM; i++, n+=2 ){
    t *= t0/i;
    sum += ( i % 4 <= 1 ? t : -t ) / n;
    if( t < fabs(sum) * zTOL ) break;
    if( i % 2 ){
      *s = sum;
      sum = *c;
    } else{
      *c = sum;
      sum = *s;
    }
  }
  if( i > Z_MAX_ITER_NUM ){
    ZITERWARN( Z_MAX_ITER_NUM );
    return false;
  }
  return true;
}

/* Fresnel integral by modified Lentz method */
static bool _zFresnelIntgLentz(double x, double *s, double *c)
{
  int i, n;
  double a, pix2;
  zComplex z1, b, ca, cd, cc, d, h, del, cs;
  bool ret = true;

  pix2 = zPI_2 * x * x;
  _zComplexCreate( &z1, 1.0, 0.0 );
  _zComplexCreate( &b, 1.0, -pix2*2 );
  _zComplexCreate( &cc, 1.0/Z_FRESNEL_FP_MIN, 0.0 );
  zComplexCDiv( &z1, &b, &d );
  _zComplexCopy( &d, &h );
  for( n=1, i=2; i<=Z_MAX_ITER_NUM; i++, n+=2 ){
    a = -n * (n+1);
    b.re += 4.0;
    _zComplexMul( &d, a, &cd );
    _zComplexAddDRC( &cd, &b );
    zComplexCDiv( &z1, &cd, &d );
    zComplexInv( &cc, &cd );
    _zComplexMulDRC( &cd, a );
    _zComplexAdd( &cd, &b, &cc );
    _zComplexCMul( &cc, &d, &del );
    zComplexCMulDRC( &h, &del );
    del.re -= 1.0;
    if( zComplexIsTiny( &del ) ) break;
  }
  if( i > Z_MAX_ITER_NUM ){
    ZITERWARN( Z_MAX_ITER_NUM );
    ret = false;
  }
  _zComplexCreateCMul( &ca, h.re, h.im, x, -x );
  _zComplexCreateCMul( &cd, cos(pix2), sin(pix2), ca.re, ca.im );
  _zComplexCreateCMul( &cs, 0.5, 0.5, 1-cd.re, -cd.im );
  *c = cs.re;
  *s = cs.im;
  return ret;
}

/* Fresnel integral */
bool zFresnelIntgPI_2(double x, double *s, double *c)
{
  double ax;
  const double xth = 1.5; /* a heuristic threshould */
  bool ret = true;

  if( ( ax = fabs( x ) ) < sqrt( Z_FRESNEL_FP_MIN ) ){ /* avoid failure of convergence */
    *s = 0;
    *c = ax;
  } else if( ax <= xth ){ /* Gilbert formula for a small argument */
    ret = _zFresnelIntgGilbert( ax, s, c );
  } else{ /* modified Lentz method for a large argument */
    ret = _zFresnelIntgLentz( ax, s, c );
  }
  if( x < 0 ){ /* antisymmetry */
    *c = -*c;
    *s = -*s;
  }
  return ret;
}

複素数使った修正Lentz法の所がちょっと見づらいですが、ご勘弁を。

Fresnel積分の一般化

さて、上記のように既に知られている数値計算法を使い、次の計算をしたいというのがこの記事の本題です。
$$
\int_0^{x} \exp {i(\phi_{0}+\phi_{1}\sigma+\phi_{2}\sigma^{2})}\mathrm{d}\sigma=C(x)+iS(x)
$$
こんなもの何に使うの？という疑問は尤もで、クロソイド補間に使います。例えば次の牧野先生の方法とか。

基本形のスケーリング

まずは次を計算する方法から考えましょう。
$$
\int_0^{x} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma=C_{2}(x;\phi_{2})+iS_{2}(x;\phi_{2})
$$
混乱を避けるために、いったん元の基本形の媒介変数を$\sigma^{\prime}$、変数を$x^{\prime}$とそれぞれ置き直します。
$$
\int_0^{x^{\prime}} e^{i\frac{\pi}{2}\sigma^{\prime 2}}\mathrm{d}\sigma^{\prime}=c(x^{\prime})+is(x^{\prime})
$$
意味は変わっていません。$\phi_{2}$は任意の実数を取り得ますが、簡単のためここでは正である$(\phi_{2}>0)$と仮定しましょう。基本形を上記に揃えるため
$$
\phi_{2}\sigma^{2}=\frac{\pi}{2}\sigma^{\prime 2}
$$
なる変数変換を行えば、
$$
\sigma^{\prime}=\sqrt{\frac{2\phi_{2}}{\pi}}\sigma,\qquad
\mathrm{d}\sigma^{\prime}=\sqrt{\frac{2\phi_{2}}{\pi}}\mathrm{d}\sigma
$$
であり、また積分範囲は
$$
\begin{align}
\sigma^{\prime}&:0\rightarrow x^{\prime}
\
\sigma&:0\rightarrow \sqrt{\frac{2\phi_{2}}{\pi}}x^{\prime}
\end{align}
$$
と対応しますので、
$$
\sqrt{\frac{2\phi_{2}}{\pi}}\int_0^{\sqrt{\frac{2\phi_{2}}{\pi}}x^{\prime}} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma=c(x^{\prime})+is(x^{\prime})
$$
が言えます。よって、
$$
\begin{align}
C_{2}(x;\phi_{2})&=\sqrt{\frac{\pi}{2\phi_{2}}}c\left(\sqrt{\frac{2\phi_{2}}{\pi}}x\right)
\
S_{2}(x;\phi_{2})&=\sqrt{\frac{\pi}{2\phi_{2}}}s\left(\sqrt{\frac{2\phi_{2}}{\pi}}x\right)
\end{align}
$$
となります。

$\phi_{2}$が負の場合は、
$$
\begin{align}
C_{2}(x;\phi_{2})=\int_0^{x}\cos\phi_{2}\sigma^{2}\mathrm{d}\sigma&=\int_0^{x}\cos\left(-\phi_{2}\sigma^{2}\right)\mathrm{d}\sigma=C_{2}(x;-\phi_{2}) \
S_{2}(x;\phi_{2})=\int_0^{x}\sin\phi_{2}\sigma^{2}\mathrm{d}\sigma&=-\int_0^{x}\sin\left(-\phi_{2}\sigma^{2}\right)\mathrm{d}\sigma=-S_{2}(x;-\phi_{2})
\end{align}
$$
とすれば良いです。

ZMでの実装は次の通りです。

/* a scaled version of Fresnel integral */
bool zFresnelIntgScale(double x, double f, double *s, double *c)
{
  bool ret;
  double r;

  if( zIsTiny( f ) ){ /* sin 0x = 0, cos 0x = 1 */
    *s = 0;
    *c = x;
    return true;
  }
  r = sqrt( zPI_2/fabs(f) );
  ret = zFresnelIntgPI_2( x/r, s, c );
  *s *= f > 0 ? r : -r;
  *c *= r;
  return ret;
}

一般化された形の平方完成

$\phi_{2}\neq 0$ならば、
$$
\begin{align}
\int_0^{x} \exp {i(\phi_{0}+\phi_{1}\sigma+\phi_{2}\sigma^{2})}\mathrm{d}\sigma
&=\int_0^{x} \exp i\left\{\phi_{2}\left(\sigma+\frac{\phi_{1}}{2\phi_{2}}\right)^{2}+\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right\}\mathrm{d}\sigma
\
&=\exp i\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right)
\int_0^{x} \exp i\phi_{2}\left(\sigma+\frac{\phi_{1}}{2\phi_{2}}\right)^{2}\mathrm{d}\sigma
\
&=\exp i\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right)
\int_{\frac{\phi_{1}}{2\phi_{2}}}^{x+\frac{\phi_{1}}{2\phi_{2}}} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma
\
&=\exp i\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right)\left(
\int_{0}^{x+\frac{\phi_{1}}{2\phi_{2}}} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma
-\int_{0}^{\frac{\phi_{1}}{2\phi_{2}}} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma
\right)
\end{align}
$$
なので、
$$
\begin{align}
\bar{c}&=\cos\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right) \
\bar{s}&=\sin\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right) \
\bar{C}(x)&=C_{2}\left(x+\frac{\phi_{1}}{2\phi_{2}};\phi_{2}\right)-C_{2}\left(\frac{\phi_{1}}{2\phi_{2}};\phi_{2}\right) \
\bar{S}(x)&=S_{2}\left(x+\frac{\phi_{1}}{2\phi_{2}};\phi_{2}\right)-S_{2}\left(\frac{\phi_{1}}{2\phi_{2}};\phi_{2}\right)
\end{align}
$$
とおけば、
$$
C(x)+iS(x)=(\bar{c}+i\bar{s})(\bar{C}(x)+i\bar{S}(x))
$$
すなわち
$$
\begin{align}
C(x)&=\bar{c}\bar{C}(x)-\bar{s}\bar{S}(x) \
S(x)&=\bar{c}\bar{S}(x)+\bar{s}\bar{C}(x)
\end{align}
$$
となります。

例外処理

$\phi_{2}=0$のときは、
$\phi_{1}\neq 0$なら
$$
\int_0^{x} \exp {i(\phi_{0}+\phi_{1}\sigma)}\mathrm{d}\sigma
=\frac{e^{i(\phi_{0}+\phi_{1}x)}-e^{i\phi_{0}}}{i\phi_{1}}
$$
より
$$
\begin{align}
C(x)&=\frac{\sin(\phi_{0}+\phi_{1}x)-\sin\phi_{0}}{\phi_{1}} \
S(x)&=-\frac{\cos(\phi_{0}+\phi_{1}x)-\cos\phi_{0}}{\phi_{1}}
\end{align}
$$
$\phi_{1}=0$なら
$$
\int_0^{x} \exp {i\phi_{0}}\mathrm{d}\sigma=xe^{i\phi_{0}}
$$
より
$$
\begin{align}
C(x)&=x\cos\phi_{0} \
S(x)&=x\sin\phi_{0}
\end{align}
$$
となります。これで全ての場合の一般化Fresnel積分が求まりました。

ZMでの実装を最後に示します。

/* generalized Fresnel integral */
bool zFresnelIntgGen(double x, double f0, double f1, double f2, double *s, double *c)
{
  double s1, c1, s2, c2;
  bool ret1 = true, ret2 = true;

  if( zIsTiny( f2 ) ){
    zSinCos( f0, &s1, &c1 );
    if( zIsTiny( f1 ) ){
      *s = x * s1;
      *c = x * c1;
    } else{
      zSinCos( f0 + f1 * x, &s2, &c2 );
      *s =-( c2 - c1 ) / f1;
      *c = ( s2 - s1 ) / f1;
    }
  } else{
    ret1 = zFresnelIntgScale( x+0.5*f1/f2, f2, &s1, &c1 );
    ret2 = zFresnelIntgScale(   0.5*f1/f2, f2, &s2, &c2 );
    s1 -= s2;
    c1 -= c2;
    zSinCos( f0 - 0.25 * f1*f1 / f2, &s2, &c2 );
    *c = c1*c2 - s1*s2;
    *s = c1*s2 + s1*c2;
  }
  return ret1 && ret2;
}

クロソイド曲線

せっかくなのでクロソイド曲線を描いてみました。(C(x),S(x))の軌跡をプロットしたものです。

$\phi_{0}$、$\phi_{1}$、$\phi_{2}$を色々変えることで中心位置、傾き、長さが変わります。