home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e031 / 3.ddi / MATHZIP2 / STARTUP / SININTEG.M < prev    next >
Encoding:
Text File  |  1991-09-23  |  5.7 KB  |  186 lines
  1.  
  2. Begin["System`"]
  3.  
  4. Unprotect[SinIntegral, CosIntegral, SinhIntegral, CoshIntegral,
  5.     FresnelC, FresnelS];
  6.  
  7. SinIntegral::usage = "SinIntegral[x] gives the sine integral
  8. Integrate[Sin[t]/t, {t, 0, x}]."
  9.  
  10. CosIntegral::usage = "CosIntegral[x] gives the cosine integral
  11. EulerGamma + Log[x] + Integrate[(Cos[t] - 1)/t, {t, 0, x}]."
  12.  
  13. SinhIntegral::usage = "SinhIntegral[x] gives the hyperbolic sine integral
  14. Integrate[Sinh[t]/t, {t, 0, x}]."
  15.  
  16. CoshIntegral::usage = "CoshIntegral[x] gives the hyperbolic cosine integral
  17. EulerGamma + Log[x] + Integrate[(Cosh[t] - 1)/t, {t, 0, x}]."
  18.  
  19. FresnelS::usage = "FresnelS[x] gives the Fresnel integral
  20. S[x] = Integrate[Sin[Pi t^2/2], {t, 0, x}]."
  21.  
  22. FresnelC::usage = "FresnelC[x] gives the Fresnel integral
  23. C[x] = Integrate[Cos[Pi t^2/2], {t, 0, x}]."
  24.  
  25. Begin["`Private`"]
  26. SinIntegral[z_]:=    (* evaluated for inexact z. *)
  27.     Module[{x, y},
  28.     x = Re[z];
  29.     y = Im[z];
  30.     Which[
  31.         z == 0,    z,
  32.         x + y < 0,    -SinIntegral[-z],
  33.         y == 0,    Im[ExpIntegralE[1,I x]] + N[Pi/2,Precision[z]] + I y/x,
  34.         x >= y,    -I/2 (ExpIntegralE[1, I z] - ExpIntegralE[1, -I z]) +
  35.                 N[Pi/2, Precision[z]],
  36.         True,    I/2 (ExpIntegralEi[-I z] + ExpIntegralE[1,-I z])
  37.     ]
  38.     ] /; (NumberQ[z] && Precision[z] < Infinity)
  39.  
  40. CosIntegral[z_]:=    (* evaluated for inexact z. *)
  41.     Module[{x, y},
  42.     x = Re[z];
  43.     y = Im[z];
  44.     Which[
  45.         z == 0,    -Infinity,
  46.         x + y < 0,    CosIntegral[-z] + If[y < 0, -I, I] N[Pi, Precision[z]],
  47.         y == 0,    -Re[ExpIntegralE[1, I x]] + I y/x,
  48.         x >= y,    -(ExpIntegralE[1, I z] + ExpIntegralE[1, -I z])/2,
  49.         True,    (ExpIntegralEi[-I z] - ExpIntegralE[1, -I z])/2 +
  50.                 I N[Pi/2, Precision[z]]
  51.     ]
  52.     ] /; (NumberQ[z] && Precision[z] < Infinity)
  53.  
  54. SinhIntegral[z_]:=    (* evaluated for inexact z. *)
  55.     Module[{shi = -I SinIntegral[I z]},
  56.     If[Head[z] === Real, Re[shi], shi]
  57.     ] /; (NumberQ[z] && Precision[z] < Infinity)
  58.  
  59. CoshIntegral[z_]:=    (* evaluated for inexact z. *)
  60.     Module[{chi},
  61.     If[z == 0, Return[-Infinity]];
  62.     chi = CosIntegral[I z] + If[Re[z] < 0 && Im[z] >= 0, 3I, -I] *
  63.         N[Pi/2, Precision[z]];
  64.     If[Head[z] === Real && z > 0, Re[chi], chi]
  65.     ] /; (NumberQ[z] && Precision[z] < Infinity)
  66.  
  67. FresnelS[z_] :=
  68.     Module[{ipiz2 = N[Pi, Precision[z]] I z^2/2},
  69.     If[Head[z] === Real,
  70.         z Im[Hypergeometric1F1[1/2, 3/2, ipiz2]],
  71.         -I z (Hypergeometric1F1[1/2, 3/2, ipiz2] -
  72.         Hypergeometric1F1[1/2, 3/2, -ipiz2])/2
  73.     ]
  74.     ] /; (NumberQ[z] && Precision[z] < Infinity)
  75.  
  76. FresnelC[z_] :=
  77.     Module[{ipiz2 = N[Pi, Precision[z]] I z^2/2},
  78.     If[Head[z] === Real,
  79.         z Re[Hypergeometric1F1[1/2, 3/2, ipiz2]],
  80.         z (Hypergeometric1F1[1/2, 3/2, ipiz2] +
  81.         Hypergeometric1F1[1/2, 3/2, -ipiz2])/2
  82.     ]
  83.     ] /; (NumberQ[z] && Precision[z] < Infinity)
  84.  
  85. SinIntegral[0]:=0
  86. CosIntegral[0]:=-Infinity
  87. SinhIntegral[0]:=0
  88. CoshIntegral[0]:=-Infinity
  89. FresnelS[0]:=0
  90. FresnelC[0]:=0
  91.  
  92. SinIntegral[DirectedInfinity[d_]] :=
  93.     Switch[d,
  94.         1,     Pi/2,
  95.         I,    DirectedInfinity[I],
  96.         -1,    -Pi/2,
  97.         -I,    DirectedInfinity[-I]
  98.     ] /; MemberQ[{1, I, -1, -I}, d]
  99.  
  100. CosIntegral[DirectedInfinity[d_]] :=
  101.     Switch[d,
  102.         1,     0,
  103.         I,    DirectedInfinity[1] + I Pi/2,
  104.         -1,    -I Pi,
  105.         -I,    DirectedInfinity[1] - I Pi/2
  106.     ] /; MemberQ[{1, I, -1, -I}, d]
  107.  
  108. SinhIntegral[DirectedInfinity[d_]] :=
  109.     Switch[d,
  110.         1,     DirectedInfinity[1],
  111.         I,    I Pi/2,
  112.         -1,    DirectedInfinity[-1],
  113.         -I,    -I Pi/2
  114.     ] /; MemberQ[{1, I, -1, -I}, d]
  115.  
  116. CoshIntegral[DirectedInfinity[d_]] :=
  117.     Switch[d,
  118.         1,     DirectedInfinity[1],
  119.         I,    I Pi/2,
  120.         -1,    DirectedInfinity[1] + I Pi,
  121.         -I,    -I Pi/2
  122.     ] /; MemberQ[{1, I, -1, -I}, d]
  123.  
  124. FresnelS[DirectedInfinity[d_]] :=
  125.     Switch[d,
  126.         1,     1/2,
  127.         I,    -I/2,
  128.         -1,    -1/2,
  129.         -I,    I/2
  130.     ] /; MemberQ[{1, I, -1, -I}, d]
  131.  
  132. FresnelC[DirectedInfinity[d_]] :=
  133.     Switch[d,
  134.         1,     1/2,
  135.         I,    I/2,
  136.         -1,    -1/2,
  137.         -I,    -I/2
  138.     ] /; MemberQ[{1, I, -1, -I}, d]
  139.  
  140. SinIntegral/: SinIntegral' := Sin[#]/# &
  141. CosIntegral/: CosIntegral' := Cos[#]/# &
  142. SinhIntegral/: SinhIntegral' := Sinh[#]/# &
  143. CoshIntegral/: CoshIntegral' := Cosh[#]/# &
  144. FresnelS/: FresnelS' := Sin[Pi/2 #^2]&
  145. FresnelC/: FresnelC' := Cos[Pi/2 #^2]&
  146.  
  147. End[]
  148.  
  149. Attributes[SinIntegral] = {Listable, ReadProtected};
  150. Attributes[CosIntegral] = {Listable, ReadProtected};
  151. Attributes[SinhIntegral] = {Listable, ReadProtected};
  152. Attributes[CoshIntegral] = {Listable, ReadProtected};
  153. Attributes[FresnelS] = {Listable, ReadProtected};
  154. Attributes[FresnelC] = {Listable, ReadProtected};
  155.  
  156. Protect[SinIntegral, CosIntegral, SinhIntegral, CoshIntegral,
  157.     FresnelS, FresnelC]
  158.  
  159. End[]
  160.  
  161. (* :Tests:
  162.  
  163. d1[z_] := Chop[1 - SinIntegral[z] / NIntegrate[Sin[t]/t, {t, 0, z},
  164.         PrecisionGoal -> 13, Method -> DoubleExponential]];
  165. d2[z_] := Chop[1 - SinhIntegral[z] / NIntegrate[Sinh[t]/t, {t, 0, z},
  166.         PrecisionGoal -> 13, Method -> DoubleExponential]];
  167. d3[z_] := Chop[1 - CosIntegral[z] / (N[EulerGamma] + Log[z] +
  168.         NIntegrate[(Cos[t] - 1)/t, {t, 0, z},
  169.             PrecisionGoal -> 13, Method -> DoubleExponential])];
  170. d4[z_] := Chop[1 - CoshIntegral[z] / (N[EulerGamma] + Log[z] +
  171.         NIntegrate[(Cosh[t] - 1)/t, {t, 0, z},
  172.             PrecisionGoal -> 13, Method -> DoubleExponential])];
  173. d5[z_] := Chop[1 - FresnelS[z] / NIntegrate[Sin[Pi/2 t^2], {t, 0, z},
  174.         PrecisionGoal -> 13, Method -> DoubleExponential]];
  175. d6[z_] := Chop[1 - FresnelC[z] / NIntegrate[Cos[Pi/2 t^2], {t, 0, z},
  176.         PrecisionGoal -> 13, Method -> DoubleExponential]];
  177.  
  178. Do[zz = Random[Complex, {-20-20I, 20+20I}];
  179.     Print[{zz, d1[zz], d2[zz], d3[zz], d4[zz], d5[zz], d6[zz]}], {20}];
  180. Do[zz = Random[Real, {-20, 20}];
  181.     Print[{zz, d1[zz], d2[zz], d3[zz], d4[zz], d5[zz], d6[zz]}], {20}];
  182. Do[zz = Random[Real, {-20, 20}]I;
  183.     Print[{zz, d1[zz], d2[zz], d3[zz], d4[zz], d5[zz], d6[zz]}], {20}];
  184.  
  185. *)
  186.