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}$を色々変えることで中心位置、傾き、長さが変わります。