Hi folks,

sadly it has been a while and I think I have to learn the things from anew how to use the NR.

So my problem is that I try to sample my monte carlo method by root finding.

The procedure goes as follows:
I have a cross section for my interaction, depending on the primary energy T and secondary energy W: \sigma(T,W).

To describe the interaction via a random number I divide the cross section by the maximum possible energy of the secondary electron
\frac{\sigma(T,W)}{\sigma(T,W_{max})} = u
where u is a random number in [0,1].

So if I put in the information of a random number, I get
\frac{\sigma(T,W)}{\sigma(T,W_{max})} = u \quad \Rightarrow\quad \frac{\sigma(T,W)}{\sigma(T,W_{max})} - u = 0

Well, this is a root finding problem, where
f(x) = \frac{\sigma(T,x)}{\sigma(T,W_{max})} - u

To use the Newton Raphson Method I have to derive my function f, leading to
f'(x) = \frac{1}{\sigma(T,W_{max})} \frac{d\sigma(T,x)}{dx}

Now, the Newton Raphson Method is based on iteration via
x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}

Now, for my problem this leads to
x_{i+1} = x_i - \frac{\sigma(T,x) - u\cdot \sigma(T,W_{max})}{\frac{d\sigma(T,x)}{dx}}

Now, this should work or not? I have the analytical description of both the cross section sigma and the derivative.
I guess I have to work as in post nr. 5 (http://www.nr.com/forum/showpost.php?p=4677&postcount=5) with the same problems.

Cheers,
Quantum.

So, I was busy tu bild a functor according to my last example, given in the link (http://www.nr.com/forum/showthread.php?t=1762&highlight=quadrature) and now for the rtnewt - Method according to the Numerical Recipes, 3rd ed. on page 459:


// IonizationStructFuncdNR.h
struct Funcd{
// constructor
Funcd(Doub T, Doub W, Doub I);
// variables
Doub t,w;
Doub S;
// total cross Section
Doub operator () (const Doub x)
{
return 0.5 * (// according to the Rudd Paper
S*F1(t)*(0.5 - 1/(2*pow(w+1,2)) - 1/(2*pow(t,2)) //
+ 1/(2*pow(t-w,2)) //
+ 2*(t+2*w-1) / ( pow(t+1,2)*sqrt(w+1)*sqrt(t-w) ) //
- 2*(t-1) / ( pow(t+1,2)*sqrt(t) ) )
+
S*F2(t)*( 1 - 1/(w+1) - 1/t + 1/(t-w) //
+ 1/(t+1) * log((t-w)/(t*(w+1))) )
);
}
// derivative
Doub df(const Doub x)
{
return
S/I*F1(t)*(1/pow(1+w,3) + 1/pow(t-w,3) //
- 1/( pow(1+w,3/2) * pow(t-w,3/2) )
+
S/I*F2(t)*(1/pow(1+w,2) + 1/pow(t-w,2) //
- 1/( (1+w) * (t-w) ) )
);
}
};
// constructor
Funcd::Funcd(Doub T, Doub W, Doub I)
{
t = T/I;
w = W/I;
S = 4*M_PI*pow(a0,2)*N*pow(Ry/I,2);
}


Now I have the problem, that after compiling I get the error message:

25: error: ?I? was not declared in this scope

Line 25 refers to the section about the derivative. But I wonder why it should not be declared. I declare it in the construre via Funcd(Doub T, Doub W, Doub I) - so what is the problem?

Cheers,
Quantum.