pinit.c File Reference


Detailed Description

We allocate and initialize a new Parratt struct, parse the XML file, and if needed, load the data.


Functions

ParrattinitParratt (char *filename)
 Returns a newly allocated Parratt struct.
void initData (struct _data *data)
 Loads the data from file.
int initMethod (char *type, char *algo)
 Initializes the fit method types.
void initResolution (struct _other *other)
 Loads the resolution values from file.


Function Documentation

Parratt* initParratt ( char *  filename  ) 

Initializes a newly allocated Parratt struct. The struct is loaded with the XML input by calling loadXML(). If the do_fit flag is "yes", then the data to be fit is loaded as well by calling initData().

00027 {
00028     unsigned int i, j;
00029     Parratt *par;
00030     
00031     /*---------------------------------
00032      *  Allocate a new Parratt struct.
00033      */ 
00034     par = (Parratt *)malloc(sizeof(Parratt));
00035     
00036     /*---------------------------------
00037      *  Set the misc variables to defaults.
00038      */
00039     par->s.doit = 0;
00040     par->s.qmin = 0.0;
00041     par->s.qmax = 0.3;
00042     par->s.numq = 200;
00043     par->s.numz = 200;
00044     for (j=0; j<64; j++) 
00045     {
00046         par->s.refl[j] = 0;
00047         par->s.prof[j] = 0;
00048     }
00049     sprintf(par->s.refl, "reflectivity.txt");
00050     sprintf(par->s.prof, "profile.txt");
00051     par->s.pntp = PARRATT;
00052     
00053     /*---------------------------------
00054      *  Zero the layers struct
00055      */
00056     par->l.num = 0;
00057     for (i=0; i<MAX_LAYER; i++)
00058     {
00059         for (j=0; j<64; j++)
00060         {
00061             par->l.nam[i][j] = 0;
00062         }
00063         par->l.idx[i] = 0;
00064         par->l.typ[i] = 0;
00065     }
00066     
00067     /*---------------------------------
00068      *  Zero the model struct
00069      */
00070     par->m.num  = 2;
00071     par->m.idx  = 0;
00072     for (i=0; i<MAX_PARAM; i++)
00073     {
00074         par->m.par[i] = 0.0;
00075         par->m.fix[i] = 1;
00076     }
00077     par->m.par[1] = 1.0;
00078 
00079     /*---------------------------------
00080      *  Zero the other struct
00081      */
00082     par->g.plz = POLARIZED_UP;
00083     par->g.rtp = CONSTANT;
00084     par->g.rsn = 0.0;
00085     par->g.rl = SINGLE;
00086     par->g.rsNum = 0;
00087     for (i=0; i<64; i++) par->g.rnam[i] = 0;
00088         
00089     /*---------------------------------
00090      *  Zero the data structs
00091      */
00092     for (i=0; i<64; i++) par->d1.nam[i] = 0;
00093     par->d1.idx = 0;
00094     par->d1.num = 0;
00095     par->d1.wgt = WEIGHT_NONE;
00096 
00097     for (i=0; i<64; i++) par->d2.nam[i] = 0;
00098     par->d2.idx = 0;
00099     par->d2.num = 0;
00100     par->d2.wgt = WEIGHT_NONE;
00101     
00102     /*---------------------------------
00103      *  Parse the XML file
00104      */
00105     loadXML(filename, par);
00106     if (par->l.num < 2) error("Not enough layers");
00107 
00108     /*---------------------------------*/
00109     /*
00110      *  The 'do' attribute flags whether to load data and par.
00111      */
00112     if (par->s.doit) 
00113     {   
00114         initData(&par->d1);
00115         if (par->s.pntp == POLARIZED)
00116         {
00117             initData(&par->d2);
00118         }
00119     }
00120 
00121     /*---------------------------------*/
00122     /*
00123      *  Get resolution values from file, if required.
00124      */ 
00125     if (par->g.rl==ARRAY) initResolution(&par->g);
00126     
00127     return par;
00128 }

void initData ( struct _data data  ) 

This function loads the data from file. It is expecting two or three column data. Data order is unimportant.

00131 {
00132     unsigned int  i, j;
00133     char          buff[120];
00134     unsigned int  ndat;
00135     double        dtmp;
00136     FILE         *fp = NULL;
00137     
00138     /*--------------------------------------------
00139      * Open the data file
00140      */
00141     fp = fopen(data->nam, "r");
00142     if (fp==NULL) error("Cannot open data file");
00143 
00144     /*--------------------------------------------
00145      *  Count the number of data points
00146      */ 
00147     ndat= 0;
00148     for (i=0; i<120; i++) buff[i] = 0;
00149     while(fgets(buff, 119, fp) != NULL)
00150     {
00151         j = sscanf(buff, "%lf %lf %lf", &dtmp, &dtmp, &dtmp);
00152         if ((j == 2) || (j == 3)) ndat++;
00153         for (i=0; i<120; i++) buff[i] = 0;
00154     }
00155     rewind(fp);
00156 
00157     /*--------------------------------------------
00158      *  Allocate memory for the data
00159      */ 
00160     data->num = ndat;
00161     data->xpt = (double *)calloc(data->num,sizeof(double));
00162     data->ypt = (double *)calloc(data->num,sizeof(double));
00163     data->sig = (double *)calloc(data->num,sizeof(double));
00164     
00165     /*--------------------------------------------
00166      *  Read in  the data
00167      */ 
00168     ndat= 0;
00169     for (i=0; i<120; i++) buff[i] = 0;
00170     while(fgets(buff, 119, fp) != NULL)
00171     {
00172         j = sscanf(buff, "%lf %lf %lf", &data->xpt[ndat], \
00173                                         &data->ypt[ndat], \
00174                                         &data->sig[ndat]);
00175         /* 
00176          *  If we want sig from file, but don't get it,
00177          *  what should we do? For now, set sig=1.
00178          */
00179         if ((j == 2) && (data->wgt == WEIGHT_ERR)) data->sig[ndat] = 1.0;
00180         if ((j == 2) || (j == 3)) ndat++;
00181         for (i=0; i<120; i++) buff[i] = 0;
00182     }
00183     fclose(fp);
00184     
00185     /*--------------------------------------------
00186      *  Set the sig variables
00187      */ 
00188     for (i=0; i<data->num; i++)
00189     {
00190         if      (data->wgt == WEIGHT_NONE) data->sig[i] = 1.0;
00191         else if (data->wgt == WEIGHT_R)    data->sig[i] = data->ypt[i];
00192         else if (data->wgt == WEIGHT_SQRT) data->sig[i] = sqrt(data->ypt[i]);
00193     }
00194     
00195     return;
00196 }

int initMethod ( char *  type,
char *  algo 
)

Default fitting type is LEAST_SQUARES.

00199 {
00200     int i = 0, j = 0;
00201     /*
00202      *  In case no functions were in the XML file,
00203      *  set some defaults.
00204      */
00205     if ((type==NULL) || (algo==NULL))
00206     {
00207         fit_method.solv_type = &lstsq_main;
00208         fit_method.lstsq_type = gsl_multifit_fdfsolver_lmder;
00209         return GSL_SUCCESS;
00210     }
00211     /*
00212      *  The solver and method types are declared here so that
00213      *  the user could choose
00214      *  "SORT OF" CHECKED FOR VALID ENTRIES
00215      */
00216     if      (strncmp(type,"MULTIDIM_MINIMIZAION", 20)==0)
00217     {
00218         fit_method.solv_type = &mdmin_main; 
00219         i = 1;
00220     }       
00221     else if (strncmp(type,"LEAST_SQUARES", 13)==0)
00222     {
00223         fit_method.solv_type = &lstsq_main;
00224         i = 2;
00225     }   
00226     else if (strncmp(type,"SIMULATED_ANNEAL", 16)==0)
00227     {
00228         fit_method.solv_type = &siman_main;
00229         i = 3;
00230     }   
00231     /* 
00232      *  Should we error out in this case?
00233      */
00234     else 
00235     {
00236         fit_method.solv_type = &lstsq_main;
00237         i = 2;
00238     }
00239 
00240     /*------------------------------------------------------*/
00241     if      (strncmp(algo,"CONJUGATE_FR", 12)==0)
00242     {
00243         fit_method.mdmin_type = gsl_multimin_fdfminimizer_conjugate_fr;
00244         j = 1;
00245     }
00246     else if (strncmp(algo,"CONJUGATE_PR", 12)==0)
00247     {
00248         fit_method.mdmin_type = gsl_multimin_fdfminimizer_conjugate_pr;
00249         j = 1;
00250     }
00251     else if (strncmp(algo,"VECTOR_BFGS", 11)==0)
00252     {
00253         fit_method.mdmin_type = gsl_multimin_fdfminimizer_vector_bfgs;
00254         j = 1;
00255     }
00256     else if (strncmp(algo,"STEEP_DESC", 10)==0)
00257     {
00258         fit_method.mdmin_type = gsl_multimin_fdfminimizer_steepest_descent;
00259         j = 1;
00260     }
00261     else if (strncmp(algo,"LMSDR", 5)==0)
00262     {
00263         fit_method.lstsq_type = gsl_multifit_fdfsolver_lmsder;
00264         j = 2;
00265     }
00266     else if (strncmp(algo,"LMDR", 4)==0)
00267     {
00268         fit_method.lstsq_type = gsl_multifit_fdfsolver_lmder;
00269         j = 2;
00270     }
00271     /* 
00272      *  Should we error out in this case?
00273      */
00274     else
00275     {
00276         fit_method.lstsq_type = gsl_multifit_fdfsolver_lmder;
00277         j = 2;
00278     }
00279 
00280     /*------------------------------------------------------*/
00281     if (((i==1) && (j==2)) || ((i==2) && (j==1)))
00282     {
00283         GSL_ERROR("Error: wrong method/algorithm pair.", GSL_EINVAL); 
00284         return GSL_EINVAL; 
00285     }
00286 
00287     return GSL_SUCCESS;
00288 }

void initResolution ( struct _other other  ) 

This function loads the resolution data from the specified file. The function expects the data in two column format: the first is the q's, the second column the associated resolutions. The data MUST be in ascending order of q for the loading to work properly. However, not all q must be specified. The listing can be thought of as defining 'q-ranges,' with each number pair defining the resolution starting at that q, and ending at the next specified q.

00291 {
00292     unsigned int i, j;
00293     char         buff[120];
00294     unsigned int ndat;
00295     double       dtmp;
00296     FILE    *fp = NULL;
00297     
00298     /*--------------------------------------------
00299      * Open the resolution file
00300      */
00301     fp = fopen(other->rnam, "r");
00302     if (fp==NULL) error("Cannot open resolution file");
00303 
00304     /*--------------------------------------------
00305      *  Count the number of data points
00306      */ 
00307     ndat= 0;
00308     for (i=0; i<120; i++) buff[i] = 0;
00309     while(fgets(buff, 119, fp) != NULL)
00310     {
00311         j = sscanf(buff, "%lf %lf %lf", &dtmp, &dtmp, &dtmp);
00312         if ((j == 2) || (j == 3)) ndat++;
00313         for (i=0; i<120; i++) buff[i] = 0;
00314     }
00315     rewind(fp);
00316     
00317     /*--------------------------------------------
00318      *  Allocate memory for the data
00319      */ 
00320     other->rsNum = ndat;
00321     other->rsArrayq   = (double *)calloc(other->rsNum,sizeof(double));
00322     other->rsArrayval = (double *)calloc(other->rsNum,sizeof(double));
00323     
00324     /*--------------------------------------------
00325      *  Read in  the data
00326      */ 
00327     ndat= 0;
00328     for (i=0; i<120; i++) buff[i] = 0;
00329     while(fgets(buff, 119, fp) != NULL)
00330     {
00331         j = sscanf(buff, "%lf %lf", &other->rsArrayq[ndat], \
00332                                         &other->rsArrayval[ndat] );
00333         if ((j == 2) || (j == 3)) ndat++;
00334         for (i=0; i<120; i++) buff[i] = 0;
00335     }
00336     fclose(fp);
00337     
00338     
00339 }


Generated on Wed Mar 14 13:24:56 2007 for Parratt by  doxygen 1.4.7