Message Boards Message Boards

0
|
4444 Views
|
4 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Ask a question about function calculation in Mathematica, thanks :)

Posted 9 years ago

Hi All, I have a question which has bothered me for a whole week. I defined some global functions and then calculated two specific function values based on them. I can plot the first one before running the calculation for the second function, but not after I running the calculation for the second function. To make my confusion clear, I copied my code below:

(*Define public parameters*)
c = 2.99792458*10^8;
kB = 1.38*10^-23;
\[Mu]0 = 4 \[Pi]*10^-7;
n1 = 1.4469;
n2 = 1.0;
a = 200*10^-9;
\[Omega][\[Lambda]_] := 2*\[Pi]*c/\[Lambda];
k0[\[Lambda]_] := 2*\[Pi]/\[Lambda];
JD[l_, x_] := BesselJ[l - 1, x] - l*BesselJ[l, x]/x;
KD[l_, x_] := -(BesselK[l - 1, x] + BesselK[l + 1, x])/2;
s11[\[Lambda]_] := (1/(h11[\[Lambda]]*a)^2 + 
     1/(q11[\[Lambda]]*a)^2)/(JD[1, h11[\[Lambda]]*a]/h11[\[Lambda]]/
       a/BesselJ[1, h11[\[Lambda]]*a] + 
     KD[1, q11[\[Lambda]]*a]/q11[\[Lambda]]/a/
      BesselK[1, q11[\[Lambda]]*a]);
f[\[Lambda]_] := (1 + s11[\[Lambda]])^2/(1 - s11[\[Lambda]])^2;
fp[\[Lambda]_] := 2*(1 + s11[\[Lambda]])/(1 - s11[\[Lambda]]);
r[x_, y_] := Sqrt[x^2 + y^2];
e = 1.602176*10^-19;
me = 9.109381*10^-31;
\[Omega]\[Alpha] = {2 \[Pi]*384.2304844685*10^12, 
   2 \[Pi]*377.107463380*10^12};(*D2 and D1 lines*)
f\[Alpha] = {0.69577, 0.34231};(*Oscillator damping coefficients*)
\[Gamma]\[Alpha] = {38.117*10^6, 36.129*10^6};(*natual line widths*)
\[Alpha]1[\[Lambda]_] := 
  e^2/me*Sum[
    f\[Alpha][[
      i]]*(\[Omega]\[Alpha][[
         i]]^2 - \[Omega][\[Lambda]]^2)/((\[Omega]\[Alpha][[
           i]]^2 - \[Omega][\[Lambda]]^2)^2 + \[Gamma]\[Alpha][[
          i]]^2*\[Omega][\[Lambda]]^2), {i, 1, 2} ];
(*Blue-detuned parameters*)
P1 = 0.029;
\[Lambda]1 = 700*10^-9;
\[Beta]11 = 10.1477*10^6;(*propagation constant*)
q11[\[Lambda]_] := Sqrt[\[Beta]11^2 - k0[\[Lambda]]^2*n2^2];
h11[\[Lambda]_] := Sqrt[k0[\[Lambda]]^2*n1^2 - \[Beta]11^2];
Din1[\[Lambda]_] := (1 - 
      s11[\[Lambda]])*(1 + (1 - s11[\[Lambda]])*\[Beta]11^2/
        h11[\[Lambda]]^2)*((BesselJ[0, 
        h11[\[Lambda]]*a])^2 + (BesselJ[1, 
        h11[\[Lambda]]*a])^2) + (1 + 
      s11[\[Lambda]])*(1 + (1 + s11[\[Lambda]])*\[Beta]11^2/
        h11[\[Lambda]]^2)*((BesselJ[2, h11[\[Lambda]]*a])^2 - 
      BesselJ[1, h11[\[Lambda]]*a]*BesselJ[3, h11[\[Lambda]]*a]);
Dout1[\[Lambda]_] := (BesselJ[1, 
     h11[\[Lambda]]*
      a])^2*((1 - 
         s11[\[Lambda]])*(1 - (1 - s11[\[Lambda]])*\[Beta]11^2/
           q11[\[Lambda]]^2)*((BesselK[0, 
           q11[\[Lambda]]*a])^2 - (BesselK[1, 
           q11[\[Lambda]]*a])^2) + (1 + 
         s11[\[Lambda]])*(1 - (1 + s11[\[Lambda]])*\[Beta]11^2/
           q11[\[Lambda]]^2)*((BesselK[2, q11[\[Lambda]]*a])^2 - 
         BesselK[1, q11[\[Lambda]]*a]*
          BesselK[3, q11[\[Lambda]]*a]))/(BesselK[1, 
      q11[\[Lambda]]*a])^2;
A1[\[Lambda]_] := 
  Sqrt[4*\[Mu]0*\[Omega][\[Lambda]]*
     P1/\[Pi]/a^2/\[Beta]11]*(Din1[\[Lambda]] + Dout1[\[Lambda]])^(-1/
    2);
Alin1[\[Lambda]_] := \[Sqrt]2*A1[\[Lambda]];
u1[\[Lambda]_] := 
  2*h11[\[Lambda]]^2/\[Beta]11^2/(1 - s11[\[Lambda]])^2;
w1[\[Lambda]_] := 
  2*q11[\[Lambda]]^2/\[Beta]11^2/(1 - s11[\[Lambda]])^2;
gin1[\[Lambda]_] := (Abs[Alin1[\[Lambda]]])^2/2/u1[\[Lambda]];
gout1[\[Lambda]_] := (Abs[
      Alin1[\[Lambda]]])^2*(BesselJ[1, h11[\[Lambda]]*a])^2/2/
     w1[\[Lambda]]/(BesselK[1, q11[\[Lambda]]*a])^2;
(*circularly-polarized*)
I1in[r_] := 
  2 gin1[\[Lambda]1]*((BesselJ[0, h11[\[Lambda]1]*r])^2 + 
     u1[\[Lambda]1]*(BesselJ[1, h11[\[Lambda]1]*r])^2 + 
     f[\[Lambda]1]*(BesselJ[2, h11[\[Lambda]1]*r])^2);
I1out[r_] := 
  2 gout1[\[Lambda]1]*((BesselK[0, q11[\[Lambda]1]*r])^2 + 
     w1[\[Lambda]1]*(BesselK[1, q11[\[Lambda]1]*r])^2 + 
     f[\[Lambda]1]*(BesselK[2, q11[\[Lambda]1]*r])^2);
F2[x_, y_] := (I1in[r[x, y]]/I1out[a]) /; Sqrt[x^2 + y^2] <= a
F2[x_, y_] := (I1out[r[x, y]]/I1out[a]) /; Sqrt[x^2 + y^2] > a
Fig6 = Plot3D[F2[x, y], {x, -4 a, 4 a}, {y, -4 a, 4 a}, 
  PlotRange -> All]
U1[x_, y_] := 
  1000*(-\[Alpha]1[\[Lambda]1]*F2[x, y]/4)*2/
    kB;(*potential energy=kB*T/2,unit in mK*)
Fig7 = Plot3D[U1[x, y], {x, a, 4*a}, {y, a, 4*a}, 
  PlotRange -> Automatic]
(*Red-detuned parameters*)
P2 = 0.03;
\[Lambda]2 = 1060*10^-9;
\[Beta]11 = 6.04866*10^6;(*propagation constant*)
q11[\[Lambda]_] := Sqrt[\[Beta]11^2 - k0[\[Lambda]]^2*n2^2];
h11[\[Lambda]_] := Sqrt[k0[\[Lambda]]^2*n1^2 - \[Beta]11^2];
Din2[\[Lambda]_] := (1 - 
      s11[\[Lambda]])*(1 + (1 - s11[\[Lambda]])*\[Beta]11^2/
        h11[\[Lambda]]^2)*((BesselJ[0, 
        h11[\[Lambda]]*a])^2 + (BesselJ[1, 
        h11[\[Lambda]]*a])^2) + (1 + 
      s11[\[Lambda]])*(1 + (1 + s11[\[Lambda]])*\[Beta]11^2/
        h11[\[Lambda]]^2)*((BesselJ[2, h11[\[Lambda]]*a])^2 - 
      BesselJ[1, h11[\[Lambda]]*a]*BesselJ[3, h11[\[Lambda]]*a]);
Dout2[\[Lambda]_] := (BesselJ[1, 
     h11[\[Lambda]]*
      a])^2*((1 - 
         s11[\[Lambda]])*(1 - (1 - s11[\[Lambda]])*\[Beta]11^2/
           q11[\[Lambda]]^2)*((BesselK[0, 
           q11[\[Lambda]]*a])^2 - (BesselK[1, 
           q11[\[Lambda]]*a])^2) + (1 + 
         s11[\[Lambda]])*(1 - (1 + s11[\[Lambda]])*\[Beta]11^2/
           q11[\[Lambda]]^2)*((BesselK[2, q11[\[Lambda]]*a])^2 - 
         BesselK[1, q11[\[Lambda]]*a]*
          BesselK[3, q11[\[Lambda]]*a]))/(BesselK[1, 
      q11[\[Lambda]]*a])^2;
A2[\[Lambda]_] := 
  Sqrt[4*\[Mu]0*\[Omega][\[Lambda]]*
     P2/\[Pi]/a^2/\[Beta]11]*(Din2[\[Lambda]] + Dout2[\[Lambda]])^(-1/
    2);
Alin2[\[Lambda]_] := \[Sqrt]2*A2[\[Lambda]];
u2[\[Lambda]_] := 
  2*h11[\[Lambda]]^2/\[Beta]11^2/(1 - s11[\[Lambda]])^2;
w2[\[Lambda]_] := 
  2*q11[\[Lambda]]^2/\[Beta]11^2/(1 - s11[\[Lambda]])^2;
gin2[\[Lambda]_] := (Abs[Alin2[\[Lambda]]])^2/2/u2[\[Lambda]];
gout2[\[Lambda]_] := (Abs[
      Alin2[\[Lambda]]])^2*(BesselJ[1, h11[\[Lambda]]*a])^2/2/
     w2[\[Lambda]]/(BesselK[1, q11[\[Lambda]]*a])^2;
I3in[r_] := 
  2 gin2[\[Lambda]2]*((BesselJ[0, h11[\[Lambda]2]*r])^2 + 
     u2[\[Lambda]2]*(BesselJ[1, h11[\[Lambda]2]*r])^2 + 
     f[\[Lambda]2]*(BesselJ[2, h11[\[Lambda]2]*r])^2);
I3out[r_] := 
  2 gout2[\[Lambda]2]*((BesselK[0, q11[\[Lambda]2]*r])^2 + 
     w2[\[Lambda]2]*(BesselK[1, q11[\[Lambda]2]*r])^2 + 
     f[\[Lambda]2]*(BesselK[2, q11[\[Lambda]2]*r])^2);
F3[x_, y_] := (I3in[r[x, y]]/I3out[a]) /; Sqrt[x^2 + y^2] <= a
F3[x_, y_] := (I3out[r[x, y]]/I3out[a]) /; Sqrt[x^2 + y^2] > a
Fig8 = Plot3D[F3[x, y], {x, -4 a, 4 a}, {y, -4 a, 4 a}, 
  PlotRange -> All]
U2[x_, y_] := 
  1000*(-\[Alpha]1[\[Lambda]2]*F3[x, y]/4)*2/
    kB;(*potential energy=kB*T/2,unit in mK*)
Fig9 = Plot3D[U2[x, y], {x, a, 4*a}, {y, a, 4*a}, 
  PlotRange -> Automatic]

I am not sure whether you need to copy the code into Mathematica to see it clearly. The problem is I can only plot out U1[x,y] before running the calculation for U2[x,y] but not after. I also attached my .nb file with this posting. Thank you very much for your help and look forward to any advice. Best, Dienter code here

Attachments:
POSTED BY: Di Chang
4 Replies

Hi All, I guess we found the problem: I should not redefine beta11. Thanks for your help. This post can be closed. Best regards, Di

POSTED BY: Di Chang

Another interesting thing is in the whole calculation, everything should be real, I do not understand how the imaginary values came out.

POSTED BY: Di Chang

After the print of Fig9

Fig9

one sees using ?U1 and FullDefinition[U1] that U1 is still well defined (seemingly). So lets test the U1 on the diagonal

In[76]:= U1[a, a]
Out[76]= 1.81085*10^-13 - 2.16846*10^-13 I

In[77]:= U1[4 a, 4 a]
Out[77]= 3.1515*10^-14 - 1.61214*10^-14 I

these small imaginary parts prevent the graphics from appearing (Graphics3D wants real values) and bring you to the idea to type

Fig10 = Plot3D[Re[U1[x, y]] + U2[x, y], {x, a, 4 a}, {y, a, 4 a}, PlotRange -> Automatic]

to see U1 - its real part actually - in the sum again:

enter image description here

It lies in your responsibility to check whether U1 had imaginary parts already right after its definition and if so, why the first picture Fig7 did print.

HTH

POSTED BY: Udo Krause

Hi Udo Krause,

Thank you very much for your reply. Yes, you are right. I checked the imaginary parts of U1 right after its definition and calculation and they are all zero. However, after U2 has been calculated, U1 starts getting imaginary part, which I do not understand why. Maybe it has something to do with they sharing the same global functions definition at the beginning? Or do you think the way I define U1 will help the computer and program keep remembering its values? I am not good at computation and I guess the problem happens because the values of U1 have not been stored there and released. Do you have any suggestions that I can modify my codes to avoid this? Again, many thanks for your help.

Best regards, Di

POSTED BY: Di Chang
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract