Hi,
Because I have been less lucky with the response to some of my postings here I'd like to state a quite simple question related to my C code for use with NLSF.
My function, parplane, has a loop that calculates an array of values that is to be used in the NLSF equation. The parameter I want to optimize is also a part of the array definition so I want the array to vary freely during the fit. Now, the array must be initialized and I have been doing that in the C code for the function. I think this makes a problem for NLSF because the array is reset each time NLSF tries to improve the fit.
My question(s) then is:
Should I move the array initialization loop to "Initialize parameters" in NLSF? Do I need to declare the array also there then? Should the array be declared both for parameter initialization and in the function code?
My code is:
#define JMAX 50
#define xacc 1e-12
#define nr_roots 100 //number of roots to be found
The function:
double parplane(double a, double q) //function for use in NLSF
{
double dSum = 0; //cumulative output set to zero
double rtn[nr_roots],x1[nr_roots],x2[nr_roots],fx[nr_roots],df[nr_roots],dx[nr_roots];
int j;
int i;
The array initialization:
for (i=0; i<nr_roots; i++)
{
x1[i] = (i+1)*PI-1;
x2[i] = (i+1)*PI+1;
rtn[i] = 0.5*(x1[i]+x2[i]);
}
//these guess for the range in which a root exists could come from a
//bracketing and bisection algorithm
i = 0;
while (i < nr_roots)
{
for (j=1;j<JMAX; j++)
{
fx[i] = rtn[i] * tan(rtn[i]) / a; //the function with roots
df[i] = tan(rtn[i]) / a + rtn[i] * (1+tan(rtn[i]) * tan(rtn[i])) / a; //first derivative of the function
dx[i] = -fx[i]/df[i];
rtn[i] += dx[i];
printf("%e\n", dx[i]);
if ((x1[i]-rtn[i])*(rtn[i]-x2[i]) < 0.0){
printf("error, jumped outside bounds");
exit(1);}
if (fabs(dx[i]) < xacc){
////printf("found root after %d attempts, at %lf\n", j, rtn[i]);
printf("The value of rtn[%d] is %lf\n", i, rtn[i]);
break;}
if (j = JMAX - 1){
printf("error - exceeded max tries no root");
i++;
exit(1);}
}
dSum += exp(-(rtn[i]^2) / (a^2)) * 2 * (1 + sin(2*rtn[i])/(2*rtn[i]))^(-1) * ((2*PI*q*a)*sin(2*PI*q*a)*cos(rtn[i])-rtn[i]*cos(2*PI*q*a)*sin(rtn[i]))^2 / ((2*PI*q*a)^2-rtn[i]^2)^2;
i++;
}
return dSum;
}
Thank yoouu!
Espen
|
Simple code!
Posted - 11/19/2003 : 08:28:27 AM
