The Museum of HP Calculators

# Fresnel Integrals for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

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

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 doesn't 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)
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
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:

FIX 9     4   XEQ "FSC"   >>>>   S(4) = 0.420515754   X<>Y   C(4) = 0.498426033   ( in 11 seconds )
and  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.

01  LBL "FRI"
02  CLA
03  STO M
04  ABS
05  PI
06  SQRT
07  X>Y?
08  GTO 04     ( a three-byte GTO )
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        ( synthetic )
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        ( synthetic )
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        ( synthetic )
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:     Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4