Muller
sample Muller's potential
potential.c
Go to the documentation of this file.
00001 /*
00002  * This work is licensed under a Creative Commons 
00003  * Attribution-Noncommercial-Share Alike 3.0 United States License.
00004  * 
00005  *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
00006  * 
00007  * You are free:
00008  * 
00009  *    to Share - to copy, distribute, display, and perform the work
00010  *    to Remix - to make derivative works
00011  * 
00012  * Under the following conditions:
00013  * 
00014  *    Attribution. You must attribute the work in the manner specified by the
00015  *    author or licensor (but not in any way that suggests that they endorse you
00016  *    or your use of the work).
00017  * 
00018  *    Noncommercial. You may not use this work for commercial purposes.
00019  * 
00020  *    Share Alike. If you alter, transform, or build upon this work, you may
00021  *    distribute the resulting work only under the same or similar license to
00022  *    this one.
00023  * 
00024  * For any reuse or distribution, you must make clear to others the license
00025  * terms of this work. The best way to do this is by including this header.
00026  * 
00027  * Any of the above conditions can be waived if you get permission from the
00028  * copyright holder.
00029  * 
00030  * Apart from the remix rights granted under this license, nothing in this
00031  * license impairs or restricts the author's moral rights.
00032  */
00033 
00034 
00046 #include <config.h>
00047 
00048 #include <stdlib.h>
00049 #include <stdio.h>
00050 #include <inttypes.h>
00051 #include <math.h>
00052 #include <assert.h>
00053 
00054 #include <jwsc/base/error.h>
00055 #include <jwsc/matrix/matrix.h>
00056 
00057 #include "potential.h"
00058 
00059 
00060 #define  ERR_BUF_LEN  256
00061 
00062 
00064 static const float A[4] = {-200.0f, -100.0f, -170.0f, 15.0f};
00065 
00066 
00068 static const float a[4] = {-1.0f, -1.0f, -6.5f, 0.7f};
00069 
00070 
00072 static const float b[4] = {0.0f, 0.0f, 11.0f, 0.6f};
00073 
00074 
00076 static const float c[4] = {-10.0f, -10.0f, -6.5f, 0.7f};
00077 
00078 
00080 static const float X[4] = {1.0f, 0.0f, -0.5f, -1.0f};
00081 
00082 
00084 static const float Y[4] = {0.0f, 0.5f, 1.5f, 1.0f};
00085 
00086 
00088 static char err_buf[ERR_BUF_LEN] = {0};
00089 
00090 
00106 Error* mullers_potential_f(float* V_out, float x, float y)
00107 {
00108     static char err_buf[256] = {0};
00109 
00110     uint8_t i;
00111     float   V;
00112 
00113     V = 0;
00114 
00115     for (i = 0; i < 4; i++)
00116     {
00117         V += A[i] * exp(a[i] * (x - X[i]) * (x - X[i]) +
00118                         b[i] * (x - X[i]) * (y - Y[i]) +
00119                         c[i] * (y - Y[i]) * (y - Y[i]));
00120     }
00121 
00122 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN
00123     if (isinf(V) || isnan(V))
00124     {
00125         snprintf(err_buf, ERR_BUF_LEN-1, 
00126                 "Muller's potential precision error: x=%e y=%e", x, y);
00127         return JWSC_ECALC(err_buf);
00128     }
00129 #endif
00130 
00131     *V_out = V;
00132 
00133     return NULL;
00134 }
00135 
00144 Error* mullers_potential_d(double* V_out, double x, double y)
00145 {
00146     uint8_t i;
00147     double  V;
00148 
00149     V = 0;
00150 
00151     for (i = 0; i < 4; i++)
00152     {
00153         V += (double)A[i] * 
00154              exp((double)a[i] * (x - (double)X[i]) * (x - (double)X[i]) +
00155                  (double)b[i] * (x - (double)X[i]) * (y - (double)Y[i]) +
00156                  (double)c[i] * (y - (double)Y[i]) * (y - (double)Y[i]));
00157     }
00158 
00159 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN
00160     if (isinf(V) || isnan(V))
00161     {
00162         snprintf(err_buf, ERR_BUF_LEN-1, 
00163                 "Muller's potential precision error: x=%e y=%e", x, y);
00164         return JWSC_ECALC(err_buf);
00165     }
00166 #endif
00167 
00168     *V_out = V;
00169 
00170     return NULL;
00171 }
00172 
00198 Error* grad_mullers_potential_f
00199 (
00200     float* dx_V_out, 
00201     float* dy_V_out, 
00202     float  x, 
00203     float  y
00204 )
00205 {
00206     uint8_t i;
00207     float   dx_V, dy_V;
00208 
00209     dx_V = 0;
00210 
00211     for (i = 0; i < 4; i++)
00212     {
00213         dx_V += A[i] * (2.0 * a[i] * (x - X[i]) + b[i] * (y - Y[i])) *
00214                 exp(a[i] * (x - X[i]) * (x - X[i]) +
00215                     b[i] * (x - X[i]) * (y - Y[i]) +
00216                     c[i] * (y - Y[i]) * (y - Y[i]));
00217     }
00218 
00219     dy_V = 0;
00220 
00221     for (i = 0; i < 4; i++)
00222     {
00223         dy_V += A[i] * (2.0 * c[i] * (y - Y[i]) + b[i] * (x - X[i])) *
00224                 exp(a[i] * (x - X[i]) * (x - X[i]) +
00225                     b[i] * (x - X[i]) * (y - Y[i]) +
00226                     c[i] * (y - Y[i]) * (y - Y[i]));
00227     }
00228 
00229 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN
00230     if (isinf(dx_V) || isnan(dx_V) || isinf(dy_V) || isnan(dy_V))
00231     {
00232         snprintf(err_buf, ERR_BUF_LEN-1, 
00233                 "Gradient of Muller's potential precision error: x=%e y=%e", 
00234                 x, y);
00235         return JWSC_ECALC(err_buf);
00236     }
00237 #endif
00238 
00239     *dx_V_out = dx_V;
00240     *dy_V_out = dy_V;
00241 
00242     return NULL;
00243 }
00244 
00258 Error* grad_mullers_potential_d
00259 (
00260     double* dx_V_out, 
00261     double* dy_V_out,
00262     double  x, 
00263     double  y
00264 )
00265 {
00266     uint8_t i;
00267     double  dx_V, dy_V;
00268 
00269     dx_V = 0;
00270 
00271     for (i = 0; i < 4; i++)
00272     {
00273         dx_V += (double)A[i] * 
00274                 (2.0 * (double)a[i] * (x - (double)X[i]) +
00275                        (double)b[i] * (y - (double)Y[i])) *
00276                 exp((double)a[i] * (x - (double)X[i]) * (x - (double)X[i]) +
00277                     (double)b[i] * (x - (double)X[i]) * (y - (double)Y[i]) +
00278                     (double)c[i] * (y - (double)Y[i]) * (y - (double)Y[i]));
00279     }
00280 
00281     dy_V = 0;
00282 
00283     for (i = 0; i < 4; i++)
00284     {
00285         dy_V += (double)A[i] * 
00286                 (2.0 * (double)c[i] * (y - (double)Y[i]) +
00287                        (double)b[i] * (x - (double)X[i])) *
00288                 exp((double)a[i] * (x - (double)X[i]) * (x - (double)X[i]) +
00289                     (double)b[i] * (x - (double)X[i]) * (y - (double)Y[i]) +
00290                     (double)c[i] * (y - (double)Y[i]) * (y - (double)Y[i]));
00291     }
00292 
00293 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN
00294     if (isinf(dx_V) || isnan(dx_V) || isinf(dy_V) || isnan(dy_V))
00295     {
00296         snprintf(err_buf, ERR_BUF_LEN-1, 
00297                 "Gradient of Muller's potential precision error: x=%e y=%e", 
00298                 x, y);
00299         return JWSC_ECALC(err_buf);
00300     }
00301 #endif
00302 
00303     *dx_V_out = dx_V;
00304     *dy_V_out = dy_V;
00305 
00306     return NULL;
00307 }
00308 
00332 Error* grad_matrix_mullers_potential_d
00333 (
00334     Matrix_d** gg_out, 
00335     double     x, 
00336     double     y
00337 )
00338 {
00339     double  dx_V, dy_V;
00340 
00341     Error* err;
00342 
00343     create_matrix_d(gg_out, 2, 2);
00344 
00345     if ((err = grad_mullers_potential_d(&dx_V, &dy_V, x, y)))
00346     {
00347         return err;
00348     }
00349 
00350     (*gg_out)->elts[0][0] = dx_V*dx_V;
00351     (*gg_out)->elts[1][1] = dy_V*dy_V;
00352     (*gg_out)->elts[0][1] = (*gg_out)->elts[1][0] = dx_V*dy_V;
00353 
00354 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN
00355     if (isinf(dx_V*dx_V) || isnan(dx_V*dx_V) || 
00356         isinf(dy_V*dy_V) || isnan(dy_V*dy_V) ||
00357         isinf(dx_V*dy_V) || isnan(dx_V*dy_V))
00358     {
00359         snprintf(err_buf, ERR_BUF_LEN-1, 
00360               "Gradient of Muller's potential matrix precision error: x=%e y=%e", 
00361                 x, y);
00362         return JWSC_ECALC(err_buf);
00363     }
00364 #endif
00365 
00366     return NULL;
00367 }
00368 
00392 Error* hessian_mullers_potential_f
00393 (
00394     Matrix_f** H_V_out, 
00395     float      x, 
00396     float      y
00397 )
00398 {
00399     uint8_t i;
00400     float*  dx2_V;
00401     float*  dy2_V;
00402     float*  dxy_V;
00403 
00404     Matrix_f* H_V;
00405 
00406     create_zero_matrix_f(H_V_out, 2, 2);
00407     H_V = *H_V_out;
00408 
00409     dx2_V = &(H_V->elts[0][0]);
00410 
00411     for (i = 0; i < 4; i++)
00412     {
00413         *dx2_V += A[i] * 
00414                   (2.0*a[i] + (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i])) *
00415                               (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i]))) *
00416                   exp(a[i] * (x - X[i]) * (x - X[i]) +
00417                       b[i] * (x - X[i]) * (y - Y[i]) +
00418                       c[i] * (y - Y[i]) * (y - Y[i]));
00419     }
00420 
00421     dy2_V = &(H_V->elts[1][1]);
00422 
00423     for (i = 0; i < 4; i++)
00424     {
00425         *dy2_V += A[i] * 
00426                   (2.0*c[i] + (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i])) *
00427                               (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i]))) *
00428                   exp(a[i] * (x - X[i]) * (x - X[i]) +
00429                       b[i] * (x - X[i]) * (y - Y[i]) +
00430                       c[i] * (y - Y[i]) * (y - Y[i]));
00431     }
00432 
00433     dxy_V = &(H_V->elts[0][1]);
00434 
00435     for (i = 0; i < 4; i++)
00436     {
00437         *dxy_V += A[i] * 
00438                   (b[i] + (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i])) *
00439                           (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i]))) *
00440                   exp(a[i] * (x - X[i]) * (x - X[i]) +
00441                       b[i] * (x - X[i]) * (y - Y[i]) +
00442                       c[i] * (y - Y[i]) * (y - Y[i]));
00443     }
00444 
00445     H_V->elts[1][0] = H_V->elts[0][1];
00446 
00447 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN
00448     if (isinf(*dx2_V) || isnan(*dx2_V) || isinf(*dy2_V) || isnan(*dy2_V) ||
00449         isinf(*dxy_V) || isnan(*dxy_V))
00450     {
00451         snprintf(err_buf, ERR_BUF_LEN-1, 
00452                 "Hessian of Muller's potential precision error: x=%e y=%e", 
00453                 x, y);
00454         return JWSC_ECALC(err_buf);
00455     }
00456 #endif
00457 
00458     return NULL;
00459 }
00460 
00472 Error* hessian_mullers_potential_d
00473 (
00474     Matrix_d** H_V_out, 
00475     double     x, 
00476     double     y
00477 )
00478 {
00479     uint8_t i;
00480     double* dx2_V;
00481     double* dy2_V;
00482     double* dxy_V;
00483 
00484     Matrix_d* H_V;
00485 
00486     create_zero_matrix_d(H_V_out, 2, 2);
00487     H_V = *H_V_out;
00488 
00489     dx2_V = &(H_V->elts[0][0]);
00490 
00491     for (i = 0; i < 4; i++)
00492     {
00493         *dx2_V += A[i] * 
00494                   (2.0*a[i] + (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i])) *
00495                               (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i]))) *
00496                   exp(a[i] * (x - X[i]) * (x - X[i]) +
00497                       b[i] * (x - X[i]) * (y - Y[i]) +
00498                       c[i] * (y - Y[i]) * (y - Y[i]));
00499     }
00500 
00501     dy2_V = &(H_V->elts[1][1]);
00502 
00503     for (i = 0; i < 4; i++)
00504     {
00505         *dy2_V += A[i] * 
00506                   (2.0*c[i] + (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i])) *
00507                               (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i]))) *
00508                   exp(a[i] * (x - X[i]) * (x - X[i]) +
00509                       b[i] * (x - X[i]) * (y - Y[i]) +
00510                       c[i] * (y - Y[i]) * (y - Y[i]));
00511     }
00512 
00513     dxy_V = &(H_V->elts[0][1]);
00514 
00515     for (i = 0; i < 4; i++)
00516     {
00517         *dxy_V += A[i] * 
00518                   (b[i] + (2.0*c[i]*(y - Y[i]) + b[i]*(x - X[i])) *
00519                           (2.0*a[i]*(x - X[i]) + b[i]*(y - Y[i]))) *
00520                   exp(a[i] * (x - X[i]) * (x - X[i]) +
00521                       b[i] * (x - X[i]) * (y - Y[i]) +
00522                       c[i] * (y - Y[i]) * (y - Y[i]));
00523     }
00524 
00525     H_V->elts[1][0] = H_V->elts[0][1];
00526 
00527 #if defined MULLER_HAVE_ISINF && defined MULLER_HAVE_ISNAN
00528     if (isinf(*dx2_V) || isnan(*dx2_V) || isinf(*dy2_V) || isnan(*dy2_V) ||
00529         isinf(*dxy_V) || isnan(*dxy_V))
00530     {
00531         snprintf(err_buf, ERR_BUF_LEN-1, 
00532                 "Hessian of Muller's potential precision error: x=%e y=%e", 
00533                 x, y);
00534         return JWSC_ECALC(err_buf);
00535     }
00536 #endif
00537 
00538     return NULL;
00539 }
00540