JSOC/proj/limbfit/apps/expfit.c - view

Return to expfit.c CVS log

Up to [Development] / JSOC / proj / limbfit / apps

File: [Development] / JSOC / proj / limbfit / apps / expfit.c (download)
Revision: 1.1, Wed Jul 21 00:10:01 2010 UTC (12 years, 8 months ago) by arta
Branch: MAIN
CVS Tags: Ver_LATEST, Ver_9-5, Ver_9-41, Ver_9-4, Ver_9-3, Ver_9-2, Ver_9-1, Ver_9-0, Ver_8-8, Ver_8-7, Ver_8-6, Ver_8-5, Ver_8-4, Ver_8-3, Ver_8-2, Ver_8-12, Ver_8-11, Ver_8-10, Ver_8-1, Ver_8-0, Ver_7-1, Ver_7-0, Ver_6-4, Ver_6-3, Ver_6-2, Ver_6-1, Ver_6-0, Ver_5-14, Ver_5-13, Ver_5-12, Ver_5-11, Ver_5-10, HEAD

Add Isabelles limbfit code to cvs

#include "expfit.h"

/*--------------------------------------------------------------------------*/
/* expfit.c -- model functions for exponential + quadratic */
/* Marcelo Emilio (c) Jan 2009 Modified from GNU Library manual */
     
int expb_f (const gsl_vector * x, void *data, 
             gsl_vector * f)
     {
       size_t n = ((struct dataF *)data)->n;
       double *y = ((struct dataF *)data)->y;
       double *t = ((struct dataF *)data)->t;
       double *sigma = ((struct dataF *) data)->sigma;
     
       double A0 = gsl_vector_get (x, 0);
       double A1 = gsl_vector_get (x, 1);
       double A2 = gsl_vector_get (x, 2);
       double A3 = gsl_vector_get (x, 3);
       double A4 = gsl_vector_get (x, 4);
       double A5 = gsl_vector_get (x, 5);
     
       size_t i;
     
       for (i = 0; i < n; i++)
         {
           /* Model Yi = A0 * exp(-z^2/2) + A3 + A4*i +A5*i*i  */
           /* z=(i-A1)/A2                                     */
           /*double t = i;*/
           double z = (t[i]-A1)/A2;
           double Yi = A0 * exp (-z*z/2.) + A3 +A4*t[i] + A5*t[i]*t[i];
           gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
         }
     
       return GSL_SUCCESS;
     }

/*--------------------------------------------------------------------------*/

int expb_df (const gsl_vector * x, void *data, 
              gsl_matrix * J)
     {
       size_t n = ((struct dataF *)data)->n;
       double *sigma = ((struct dataF *) data)->sigma;
       double *t = ((struct dataF *)data)->t;

       double A0 = gsl_vector_get (x, 0);
       double A1 = gsl_vector_get (x, 1);
       double A2 = gsl_vector_get (x, 2);
       //double A3 = gsl_vector_get (x, 3);
       //double A4 = gsl_vector_get (x, 4);
       //double A5 = gsl_vector_get (x, 5);

       size_t i;
     
       for (i = 0; i < n; i++)
         {
           /* Jacobian matrix J(i,j) = dfi / dxj, */
           /* where fi = (Yi - yi)/sigma[i],      */
           /* Model Yi = A0 * exp(-z^2/2) + A3 + A4*i +A5*i*i  */
           /* z=(i-A1)/A2                         */
           /* and the xj are the parameters (A0,A1,A2,A3,A4,A5) */
           /*double t = i;*/
           double s = sigma[i];
           double z = (t[i]-A1)/A2;
           double e = exp(-z*z/2.);
           gsl_matrix_set (J, i, 0, e/s); 
           gsl_matrix_set (J, i, 1, A0/A2 * z * e/s);
           gsl_matrix_set (J, i, 2, A0/A2 * z * z * e/s);
           gsl_matrix_set (J, i, 3, 1/s);
           gsl_matrix_set (J, i, 4, t[i]/s);
           gsl_matrix_set (J, i, 5, t[i] * t[i]/s);

         }
       return GSL_SUCCESS;
     }
     
int expb_fdf (const gsl_vector * x, void *data,
               gsl_vector * f, gsl_matrix * J)
     {
       expb_f (x, data, f);
       expb_df (x, data, J);
     
       return GSL_SUCCESS;
     }

Karen Tian