0
$\begingroup$

I have used the following code to evaluate an integral (val) numerically

mu = 0.0173262004;(*attenuation coefficient for E = ?*) 
k = (mu - 0.00324543007)/(0.00324543007); 
h = 150.; 
sz = 50.; 
eg = 0.7; 
qx = 1.0; 
u = 1.0; 
the = (22.5 Pi)/180; 
f1 = 1/(l^2 + z^2); 
f2 = 1 + k mu (l^2 + z^2)^0.5; 
f3 = Exp[-(z - h)^2/(2 sz sz)]; 
f4 = Exp[-(z + h)^2/(2 sz sz)]; 
f5 = Exp[-mu (l^2 + z^2)]; 
val = NIntegrate[
  2 Pi l f1 (f2 (f3 + f4) f5), {l, 0, Infinity}, {z, 0, Infinity}] 
dr[r_] := 0.0404 0.00324543007 eg qx val/((2 Pi)^0.5 u r sz the) 
LogPlot[dr[r], {r, 0, 10000}]

Now, i will be using the parameter (sz) as a function of r. For example, sz=0.26*r^0.69; How to do this? Thanking you in advance

Improve this question
$\endgroup$
  • $\begingroup$ Do you mean you wish to define sz as a function of r instead of a constant (as it is now)? $\endgroup$ – honeste_vivere Sep 29 '14 at 13:43
  • $\begingroup$ Yes, exactly... $\endgroup$ – Anand Sep 30 '14 at 12:40
  • $\begingroup$ Then define it as sz[r_]:=0.26*r^0.69 and when calling it within something else, use sz[r], replacing r with whatever variable is relevant (r_ is a pattern so the input need not be explicitly the letter r). $\endgroup$ – honeste_vivere Sep 30 '14 at 13:20
1
$\begingroup$

You can define your val now as a function of r, rather than as a constant. Then you will be able to operate with this function. I will show it within a short example, and you can then implement into your code by analogy. Consider an integral

sz=5;    
int=NIntegrate[x^2*Exp[-sz*x^2], {x, 0, Infinity}]
    (*   0.0396333  *)

Now let us replace this sonstant sz by a function sz = 0.25*r^0.69:

int[r_] := NIntegrate[x^2*Exp[-0.25*r^0.69 x^2], {x, 0, Infinity}];

Now one can do with this function whatever he needs. Let us plot it, for example:

Plot[int[r], {r, 1, 5}, AxesLabel -> {r, int}]

giving this:

enter image description here Have fun!

Improve this answer
$\endgroup$
1
$\begingroup$

Corrected per input from Michael E2

mu = 0.0173262004;
(*attenuation coefficient for E=?*)
k = (mu - 0.00324543007)/(0.00324543007);
h = 150;
sz[r_] = 0.26*r^0.69;
eg = .7;
qx = 1.;
u = 1.;
the = (22.5 Pi)/180;
f1 = 1/(l^2 + z^2);
f2 = 1 + k mu (l^2 + z^2)^0.5;
f3[r_] = Exp[-(z - h)^2/(2 sz[r] sz[r])];
f4[r_] = Exp[-(z + h)^2/(2 sz[r] sz[r])];
f5 = Exp[-mu (l^2 + z^2)];

val[r_?NumericQ] := NIntegrate[
  Evaluate[2 Pi l f1 (f2 (f3[r] + f4[r]) f5)],
  {l, 0, Infinity}, {z, 0, Infinity}]

dr[r_?NumericQ] :=
 0.0404*0.00324543007 eg qx val[r]/((2 Pi)^(1/2) u r sz[r] the) //
  Evaluate

logPlot = LogPlot[dr[r], {r, 200, 10000},
  PlotRange -> All,
  PlotPoints -> 51]

enter image description here

dr[200]

8.80707*10^-45

Improve this answer
$\endgroup$
  • $\begingroup$ dr[100] gives errors for me. Note that LogPlot effectively uses Block[{r = 100},..] etc. $\endgroup$ – Michael E2 Oct 1 '14 at 1:56
  • $\begingroup$ The problem with dr[100] is that the parameter r_?NumericQ is localized and does not refer to the global variable r that appears in sz and its dependents. LogPlot uses Block to set r, which sets the global r. In dr[100], the local function parameter r is 100, but the global r remains an undefined variable. The proper fix is to make sz, f3, and f4 functions of r. There are shorter workarounds, but they're kludgy (e.g., dr[r0_?NumericQ] := Block[{r = r0}, <code>] and val := NIntegrate[...]). $\endgroup$ – Michael E2 Oct 4 '14 at 12:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.