hp41programs

Fresnel

Fresnel Integrals for the HP-41


Overview
 

 1°) Ascending Series
 2°) Asymptotic Series
 3°) Ascending Series & Complex Continued Fraction
 
 

1°) Ascending Series
 

-Fresnel integrals are defined by:    S(x) = §0x  sin(pi.t2/2).dt      and      C(x) = §0x  cos(pi.t2/2).dt

Formulae:

     S(x) = SUMn=0,1,2,.....  (-1)n (pi/2)2n+1 x4n+3 / ((2n+1)!(4n+3))

     C(x) = SUMn=0,1,2,.....  (-1)n (pi/2)2n x4n+1 / ((2n)!(4n+1))
 

Data Registers:  R00 = x ; R01 to R03: temp
Flags: /
Subroutines: /
 
 

01  LBL "SX"
02  STO 00
03  3
04  Y^X
05  PI
06  6
07  /
08  *
09  1
10  STO 02
11  GTO 00
12  LBL "CX"
13  STO 00
14  0
15  STO 02
16  LBL 00
17  RCL 00
18  X^2
19  PI
20  *
21  2
22  /
23  X^2
24  CHS
25  STO 03    
26  R^
27  STO 01
28  LBL 01
29  RCL 01
30  RCL 02
31  2
32  +
33  STO 02
34  ST+ X
35  3
36  -
37  ST* Y     
38  4
39  +
40  RCL 02
41  ST* Y
42  DSE X
43  *
44  /
45  RCL 03
46  *
47  STO 01    
48  +
49  X#Y?
50  GTO 01
51  END
 

 
  ( 69 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x      S(x) or C(x)

Example:

    2  XEQ "SX"  >>>>  S(2) = 0.343415677    ( in 14 seconds )
    2  XEQ "CX" >>>>  C(2) = 0.488253412      --------------
 

Note:   This method does not give accurate results for large arguments:

    S(3) and C(3) are correct to 5 decimals
 S(3.5) and C(3.5) -----------  4 --------
    S(4) and C(4) are not exact even to 1D!
 

2°) Asymptotic Series ( large arguments )
 

Formulae:       S(x) = 1/2 - f(x).cos(pi.x2/2) - g(x).sin(pi.x2/2)     &     C(x) = 1/2 + f(x).sin(pi.x2/2) - g(x).cos(pi.x2/2)

   where               f(x) ~  ( 1 - 1*3/(pi.x2)2 + 1*3*5*7/(pi.x2)4 - ....... ) / (pi.x)
     and                g(x) ~  ( 1 - 1*3*5/(pi.x2)2 + 1*3*5*7*9/(pi.x2)4 - ....... ) / (pi2.x3)
 

Data Registers:   R00 = x ; R01 to R04: temp
Flags: /
Subroutines: /
 
 

01  LBL "FSC"
02  STO 00
03  X^2
04  PI
05  *
06  STO 03
07  SIGN
08  STO 01
09  LASTX
10  1/X
11  STO 02
12  XEQ 01
13  STO 04
14  1
15  CHS
16  STO 01     
17  CHS
18  STO 02
19  XEQ 01
20  PI
21  RCL 00
22  *
23  ST/ 04
24  /
25  RCL 03
26  2
27  /
28  STO 03
29  X<>Y
30  RAD
31  P-R
32  RCL 03     
33  RCL 04
34  P-R
35  ST- T
36  RDN
37  +
38  .5
39  ST+ Z
40  X<>Y
41  -
42  DEG
43  RTN
44  LBL 01
45  RCL 02     
46  RCL 01
47  2
48  +
49  ST* Y
50  LASTX
51  +
52  STO 01
53  *
54  RCL 03
55  X^2
56  /
57  CHS
58  STO 02     
59  +
60  X#Y?
61  GTO 01
62  END
 

 
   ( 80 bytes / SIZE 005 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             /           C(x)
           X             x           S(x)

Examples:

     10   XEQ FSC  >>>>    S(10) = 0.468169979    X<>Y    C(10) =  0.499898695    ( in 5 seconds )
     4.1  XEQ FSC  >>>>   S(4.1) = 0.475798258    X<>Y    C(4.1) = 0.573695632    ( in 14 seconds )
 

-The series are already divergent too soon for  x = 4   but if you add     RND  X<>Y   RND   X<>Y   after line 59, you'll get:

   in    FIX 9     4   XEQ "FSC"   >>>>   S(4) = 0.420515754   X<>Y   C(4) = 0.498426033   ( in 11 seconds )
   in    FIX 5     3   XEQ "FSC"   >>>>   S(3) = 0.4963140       X<>Y   S(3) = 0.6057203   ( exact values are  0.4963130  and  0.6057208 to 7D )
 

3°) Ascending Series & Continued Fraction
 

-The following program uses the same series expansions as "SX" & "CX"  if  | x |  <  SQRT(PI)
-Otherwise, a complex continued fraction is employed:

          C(x) + i.S(x) = ((1+i)/2) erf z   with   z = (1-i).(x/2).(pi)1/2

 and   ( 1 - erf z ) exp z2 =  (pi) -1/2 ( 1/(z+0.5/(z+1/(z+1.5/(z+2/(z+ .... ))))) )
 
 

Data Registers: /
Flags: /
Subroutines: /

-Synthetic registers M N O P may of course be replaced by registers R01 R02 R03 R00  ( for instance )
-In this case, delete line 122 , replace lines 110-111 by 3 , add  CLX  STO 02  after line 93 and delete line 02.
-Line 08 is a three-byte GTO 04
 
 

  01  LBL "FRI"
  02  CLA
  03  STO M
  04  ABS
  05  PI
  06  SQRT
  07  X>Y?
  08  GTO 04
  09  DEG
  10  *
  11  24
  12  STO N
  13  CLX
  14  2
  15  /
  16  STO O
  17  ENTER^
  18  LBL 01
  19  X<>Y
  20  CHS
  21  STO Z
  22  X^2
  23  RCL Y
  24  X^2
  25  +
  26  DSE N
  27  X<0?
  28  GTO 02
  29  RCL N
  30  X<>Y
  31  ST+ X
  32  /
  33  ST* Z
  34  *
  35  RCL O
  36  ST+ Z
  37  +
  38  GTO 01    
  39  LBL 02
  40  PI
  41  SQRT
  42  *
  43  ST/ Z
  44  /
  45  R-P
  46  RCL M
  47  X^2
  48  90
  49  *
  50  R^
  51  -
  52  X<>Y
  53  CHS
  54  P-R
  55  1
  56  +
  57  STO Z
  58  X<>Y
  59  ST+ Z
  60  -
  61  2
  62  ST/ Z
  63  /
  64  GTO 05    
  65  LBL 03
  66  X<> P
  67  PI
  68  2
  69  ST+ N
  70  /
  71  RCL M
  72  X^2
  73  *
  74  X^2
  75  *
  76  RCL N
  77  ENTER^
  78  DSE X
  79  *
  80  /
  81  CHS
  82  STO P
  83  RCL N
  84  ST+ X
  85  ISG X
  86  CLX
  87  /
  88  X<>Y
  89  ST+ Y
  90  X#Y?
  91  GTO 03    
  92  RTN
  93  LBL 04
  94  X<>Y
  95  STO P
  96  ENTER^
  97  XEQ 03
  98  STO O
  99  1
100  STO N
101  RCL M
102  ABS
103  3
104  Y^X
105  PI
106  *
107  2
108  /
109  STO P
110  PI
111  INT
112  /
113  ENTER^
114  XEQ 03    
115  RCL O
116  LBL 05
117  RCL M
118  SIGN
119  ST* Y
120  ST* Z
121  RDN
122  CLA
123  END

 
   ( 178 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             /          S(x)
           X             x          C(x)
           L             /            x
 

Examples:

    1.5  XEQ "FRI"  >>>>   C(1.5) = 0.445261176      X<>Y     S(1.5) = 0.697504959         ( in 22 seconds )
      4       R/S          >>>>     C(4)  = 0.498426033      X<>Y       S(4)  = 0.420515754         ( in 19 seconds )
 

-Thus, "FRI" produces accurate results for all arguments.
-However, "FSC" runs faster if x is very large.
 
 

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4