Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

mnbrak.c

Go to the documentation of this file.
00001 
00002 /*
00003  *   THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE
00004  *   CHECKED IT OUT USING THE COMMAND CHECKOUT.
00005  *
00006  *    $Id: mnbrak_8c-source.html 2161 2006-05-19 16:55:03Z paulf $
00007  *
00008  *    Revision history:
00009  *     $Log$
00009  *     Revision 1.1  2006/05/19 16:55:02  paulf
00009  *     first inclusion
00009  *
00010  *     Revision 1.1  2000/02/14 18:51:48  lucky
00011  *     Initial revision
00012  *
00013  *
00014  */
00015 
00016 #include <math.h>
00017 
00018 #define GOLD 1.618034
00019 #define GLIMIT 100.0
00020 #define TINY 1.0e-20
00021 #define MAX(a,b) ((a) > (b) ? (a) : (b))
00022 #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
00023 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
00024 
00025 void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb, float *fc, float(*func)(float) )
00026 {
00027         float ulim,u,r,q,fu,dum;
00028 
00029         *fa=(*func)(*ax);
00030         *fb=(*func)(*bx);
00031         if (*fb > *fa) {
00032                 SHFT(dum,*ax,*bx,dum)
00033                 SHFT(dum,*fb,*fa,dum)
00034         }
00035         *cx=(float)((*bx)+GOLD*(*bx-*ax));
00036         *fc=(*func)(*cx);
00037         while (*fb > *fc) {
00038                 r=(*bx-*ax)*(*fb-*fc);
00039                 q=(*bx-*cx)*(*fb-*fa);
00040                 u=(float) ((*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
00041                         (2.0*SIGN(MAX(fabs(q-r),TINY),q-r)));
00042                 ulim=(float)((*bx)+GLIMIT*(*cx-*bx));
00043                 if ((*bx-u)*(u-*cx) > 0.0) {
00044                         fu=(*func)(u);
00045                         if (fu < *fc) {
00046                                 *ax=(*bx);
00047                                 *bx=u;
00048                                 *fa=(*fb);
00049                                 *fb=fu;
00050                                 return;
00051                         } else if (fu > *fb) {
00052                                 *cx=u;
00053                                 *fc=fu;
00054                                 return;
00055                         }
00056                         u=(float)((*cx)+GOLD*(*cx-*bx));
00057                         fu=(*func)(u);
00058                 } else if ((*cx-u)*(u-ulim) > 0.0) {
00059                         fu=(*func)(u);
00060                         if (fu < *fc) {
00061                                 SHFT(*bx,*cx,u,(float)(*cx+GOLD*(*cx-*bx)))
00062                                 SHFT(*fb,*fc,fu,(*func)(u))
00063                         }
00064                 } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
00065                         u=ulim;
00066                         fu=(*func)(u);
00067                 } else {
00068                         u=(float)((*cx)+GOLD*(*cx-*bx));
00069                         fu=(*func)(u);
00070                 }
00071                 SHFT(*ax,*bx,*cx,u)
00072                 SHFT(*fa,*fb,*fc,fu)
00073         }
00074 }
00075 
00076 #undef GOLD
00077 #undef GLIMIT
00078 #undef TINY
00079 #undef MAX
00080 #undef SIGN
00081 #undef SHFT

Generated on Tue May 6 09:16:05 2003 for Earthworm Libs by doxygen1.3-rc3