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