Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2012

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: can somebody help ?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg127339] Re: can somebody help ?
  • From: Bob Hanlon <hanlonr357 at gmail.com>
  • Date: Wed, 18 Jul 2012 01:38:32 -0400 (EDT)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • Delivered-to: mathgroup-newout@smc.vnet.net
  • Delivered-to: mathgroup-newsend@smc.vnet.net
  • References: <20120717053225.D4CB9687F@smc.vnet.net>

Presumably, you want to localize more of the variables inside the module.

Don't use Print to return a result.

j[l_, z_] := Sqrt[\[Pi]/2 z] BesselJ[l, z];
n[l_, z_] := Sqrt[\[Pi]/2 z] BesselY[l, z];

HKrScattering[l_, En_, nsteps_] :=
 Module[
  {meV, en, \[Angstrom], \[Sigma], \[Epsilon], V, \[Mu], \[Lambda],
 k, s, x, xm, xn, xp, \[Psi], \[Psi]m, \[Psi]n, \[Psi]p,
 F, Fn, Fp, Fm, scat, r1, r2, u1, u2,
 K, j1, j2, n1, n2, \[Delta]},
  meV = 1/(27.3*1000);
  en = En*meV;
  \[Angstrom] = 1/0.52;
  \[Sigma] = 3.57 \[Angstrom];
  \[Epsilon] = 5.9 meV;
  V[r_] = \[Epsilon] ((\[Sigma]/r)^12 - 2 (\[Sigma]/r)^6);
  \[Mu] = (mH mKr)/(mH + mKr) /.
    {mH -> 1836, mKr -> 1836*136};
  \[Lambda] = 2 \[Pi]/Sqrt[2 \[Mu] en];
  k = 2 \[Pi]/\[Lambda];
  (*Print["\[Lambda]= ",\[Lambda]]*);
  s = 5 \[Lambda]/nsteps;
  x = {0.5 \[Sigma], 0.5 \[Sigma] + s};
  \[Psi] = {0, 0.0001};
  F[r_] = 2 \[Mu] (V[r] + 1/(2 \[Mu] r^2) l (l + 1) - en);
  Do[\[Psi]m = \[Psi][[n - 1]];
   \[Psi]n = \[Psi][[n]];
   xm = x[[n - 1]];
   xn = x[[n]];
   xp = xn + s;
   x = Append[x, xp];
   Fn = F[xn];
   Fp = F[xp];
   Fm = F[xm];
   \[Psi]p = (2 \[Psi]n - \[Psi]m +
       s^2/12 (10 Fn \[Psi]n + Fm \[Psi]m))/(1 - s^2/12 Fp);
   \[Psi] = Append[\[Psi], \[Psi]p],
   {n, 2, nsteps}];
  \[Psi] = \[Psi]/Max[\[Psi]];
  scat = Transpose[{x, \[Psi]}];
  {r1, u1} = scat[[-1]];
  {r2, u2} = scat[[-2]];
  K = (r1 u2)/(r2 u1);
  j1 = j[l, k r1];
  j2 = j[l, k r2];
  n1 = n[l, k r1];
  n2 = n[l, k r2];
  \[Delta] = ArcTan[(K j1 - j2)/(K n1 - n2)]]

Table[{l, HKrScattering[l, 2, 1000]}, {l, 0, 5}] // Grid

0	1.47897
1	-0.356374
2	0.66522
3	1.40513
4	-1.32016
5	-1.44998



Bob Hanlon


On Tue, Jul 17, 2012 at 1:32 AM, raj kumar <rajesh7796gm at gmail.com> wrote:
> dear esteemed members of the forum,
> hello from malaysia
> i have this workable code for a scattering process called
> HKrScattering.See below.Given a particular value of say l=1 and an
> energy value of say En =2, and nsteps=1000,  the code gives the
> correct result for the  phase shift.
> however, i would like to generate a set of values for the phase shift
> for different values of l, say for l=0 to l=5 for En=2.
>
> I have tried using a do loop for this purpose but have not been able
> to get correct results so far.
>
> can anybody please help?
> will appreciate your contribution enormously!
> cheers
> rk
>
>
> j[l_, z_] := Sqrt[\[Pi]/2 z] BesselJ[l, z];
> n[l_, z_] := Sqrt[\[Pi]/2 z] BesselY[l, z];
>
> HKrScattering[l_, En_, nsteps_] :=
>  Module[{meV, \[Angstrom], \[Sigma], \[Mu], \[Lambda], k, s,
>    scat, \[Delta]},
>
>   meV =  1./(27.3 * 1000.);
>   en = En*meV;
>   \[Angstrom] =  1/0.52 ;
>   \[Sigma] = 3.57 \[Angstrom];
>   \[Epsilon] = 5.9 meV ;
>   V[r_] := \[Epsilon] ( (\[Sigma]/r)^12 - 2 (\[Sigma]/r)^6);
>   \[Mu] = (mH  mKr)/(mH + mKr) /. {mH -> 1836 , mKr -> 1836*136};
>   \[Lambda] = 2 \[Pi]/Sqrt[2 \[Mu] en] // N;
>   k = 2 \[Pi]/\[Lambda];
>   (*Print["\[Lambda]= ",\[Lambda]]*);
>   s = 5 \[Lambda]/nsteps;
>   x = {0.5 \[Sigma], 0.5 \[Sigma] + s};
>   \[Psi] = {0.0, 0.0001};
>   F[r_] := 2 \[Mu] ( V[r] + 1/(2 \[Mu] r^2) l (l + 1) - en);
>   Do[
>    \[Psi]m = \[Psi][[n - 1]];
>    \[Psi]n = \[Psi][[n]];
>    xm = x[[n - 1]];
>    xn = x[[n]];
>    xp = xn + s;
>    x = Append[x, xp];
>    Fn = F[xn];
>    Fp = F[xp];
>    Fm = F[
>      xm]; \[Psi]p  = (
>       2. \[Psi]n - \[Psi]m +
>        s^2/12. (10. Fn  \[Psi]n + Fm  \[Psi]m))/(1. - s^2/12. Fp);
>    \[Psi] = Append[\[Psi], \[Psi]p],
>    {n , 2, nsteps}
>    ] ;
>   \[Psi] = \[Psi]/Max[\[Psi]];
>   scat = Transpose[{x, \[Psi]}];
>   nmax = Length[scat];
>   {r1, u1} = scat[[nmax]];
>   {r2, u2} = scat[[nmax - 1]];
>   K = (r1 u2)/(r2 u1);
>
>   j1 = j[l, k r1];
>   j2 = j[l, k r2];
>   n1 = n[l, k r1];
>   n2 = n[l, k r2];
>   \[Delta] =
>    ArcTan[ (K j1 - j2)/(K n1 - n2)];(*finding the phase shift*)
>   Print[\[Delta]]
>   (*Return[{scat,\[Delta]}*)]
>



-- 
Bob Hanlon



  • Prev by Date: Evaluating complicated integral numerically
  • Next by Date: Re: Differential Equation: Not getting result
  • Previous by thread: can somebody help ?
  • Next by thread: Plotting: how does one reverse the time axis tick mark values?