00001 #ifndef NR_MNBRAK_H
00002 #define NR_MNBRAK_H
00003
00004 #include <math.h>
00005 #include "utils.h"
00006
00007 namespace nr {
00008
00009
00016 template <class Real, class Func>
00017 void mnbrak
00018 (
00019 Real *ax, Real *bx, Real *cx,
00020 Real *fa, Real *fb, Real *fc,
00021 Func* func
00022 )
00023 {
00024 const Real GOLD = 1.618034;
00025 const Real GLIMIT = 100.0;
00026 const Real TINY = 1.0e-20;
00027 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
00028 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
00029
00030
00031 Real ulim,u,r,q,fu,dum;
00032
00033 *fa=(*func)(*ax);
00034 *fb=(*func)(*bx);
00035 if (*fb > *fa) {
00036 SHFT(dum,*ax,*bx,dum)
00037 SHFT(dum,*fb,*fa,dum)
00038 }
00039 *cx=(*bx)+GOLD*(*bx-*ax);
00040 *fc=(*func)(*cx);
00041 while (*fb > *fc) {
00042 r=(*bx-*ax)*(*fb-*fc);
00043 q=(*bx-*cx)*(*fb-*fa);
00044 u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/(2.0*SIGN(MAX(fabs(q-r),TINY),q-r));
00045 ulim=(*bx)+GLIMIT*(*cx-*bx);
00046 if ((*bx-u)*(u-*cx) > 0.0) {
00047 fu=(*func)(u);
00048 if (fu < *fc) {
00049 *ax=(*bx);
00050 *bx=u;
00051 *fa=(*fb);
00052 *fb=fu;
00053 return;
00054 } else if (fu > *fb) {
00055 *cx=u;
00056 *fc=fu;
00057 return;
00058 }
00059 u=(*cx)+GOLD*(*cx-*bx);
00060 fu=(*func)(u);
00061 } else if ((*cx-u)*(u-ulim) > 0.0) {
00062 fu=(*func)(u);
00063 if (fu < *fc) {
00064 SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
00065 SHFT(*fb,*fc,fu,(*func)(u))
00066 }
00067 } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
00068 u=ulim;
00069 fu=(*func)(u);
00070 } else {
00071 u=(*cx)+GOLD*(*cx-*bx);
00072 fu=(*func)(u);
00073 }
00074 SHFT(*ax,*bx,*cx,u)
00075 SHFT(*fa,*fb,*fc,fu)
00076 }
00077 #undef SHFT
00078 #undef SIGN
00079 }
00080
00081
00082 }
00083
00084 #endif
00085
Documentation 
1.3.6