HP 20S and HP 21S Great Circle

HP 20S and HP 21S Great Circle

The following program calculates the distance between two places on Earth (or any other planet) given the coordinates latitude (λ, east is positive) and longitude (ϕ, north is positive).

Inputs:

R1:  longitude 1, R2:  latitude 1

R4:  longitude 2, R5:  latitude 2

Enter the coordinates in DD.MMSSSS.

R1, R4: ϕ1, ϕ2;   R2, R5:  λ1, λ2 

Distance, in miles, is stored in R0.  Degrees mode is set.

Formula:

distance = acos( sin ϕ1 * sin ϕ2 + cos ϕ1 * cos ϕ2 * cos (λ2 – λ1) )* 3958.75 * π/180

On the HP 20S and HP 21S, can multiply by π/180 by executing the >RAD function.

HP 20S and HP 21S Program:  Great Circle

STEP	CODE	KEY
01	61, 41, A	LBL A
02	61, 23	DEG
03	22, 1	RCL 1
04	51, 54	>HR
05	21, 1	STO 1
06	23	SIN
07	55	*
08	22, 4	RCL 4
09	51, 54	>HR
10	21, 4	STO 4
11	23	SIN
12	75	+
13	22, 1	RCL 1
14	24	COS
15	55	*
16	22, 4	RCL 4
17	24	COS
18	55	*
19	33	(
20	22, 2	RCL 2
21	51, 54	>HR
22	21, 2	STO 2
23	65	-
24	22, 5	RCL 5
25	51, 54	>HR
26	21, 5	STO 5
27	34	)
28	24	COS
29	74	=
30	51, 24	ACOS
31	55	*
32	3	3
33	9	9
34	5	5
35	8	8
36	73	.
37	7	7
38	5	5
39	74	=
40	61, 55	>RAD
41	21, 0	STO 0
42	61, 26	RTN

Example 1:

Los Angeles to Rome:

Los Angeles (ϕ1 = 34°13’, λ1 = -118°15’)

Rome (ϕ2 = 41°15’, λ2 = 12°30’)

Result:  6322.2196 mi

Example 2:

Dublin to Las Vegas:

Dublin (ϕ1 = 53°20’52”, λ1 = -6°15’35”)

Las Vegas (ϕ2 = 36°10’30”, λ2 = -115°08’11”)

Result:  4938.7520 mi

Eddie

This blog is property of Edward Shore, 2017.

Fun With the HP 71B III

Fun With the HP 71B III

Links to other HP 71B Programs:

Fun with the 71B:
http://edspi31415.blogspot.com/2012/06/fun-with-hp-71b.html

Fun with the 71B II:
http://edspi31415.blogspot.com/2016/06/fun-with-hp-71b-ii.html

71B: Cubic Polynomials:
http://edspi31415.blogspot.com/2012/06/cubic-formula-basic-program-hp-71b.html

3x3 Matrices: Determinant, Inverse, 3x3 Linear Systems
EWS 6/29/2016

The program MATX3 calculates:

1. The determinant and (if possible), the inverse of a 3x3 matrix M.
2. The solution to a 3x3 linear system: Mq=D. The determinant of M will also be displayed.

If det(M) = 0, then the matrix is singular and execution stops.

The matrix M is broken into three columns (3x1 arrays): [ M ] = [ A | B | C ].

Hence M = [[ A1 B1 C1 ] [ A2 B2 C2 ] [ A3 B3 C3 ]]

Other variables used:
E = det(M)
I = M^-1. Unlike M, I will be a 3 x 3 array.
R, K, S, H: other variables used

Program MATX3 (767 bytes)
10 DESTROY A,B,I,C,R,K,S,H,D,Q
11 DISP “1. DET/INV  2. 3x3” @ WAIT 2
12 INPUT “1. D/I 2. SYS:”; H
13 DIM A(3),B(3),C(3),I(3,3),D(3)
14 OPTION BASE 1
20 FOR K=1 TO 3
21 DISP “ROW “;K @ WAIT 1
22 INPUT “A:”; A(K)
24 INPUT “B:”; B(K)
26 INPUT “C:”; C(K)
28 IF H=2 THEN INPUT “D:”; D(K)
30 NEXT K
40 DEF FND(X,Y,Z,T)=X*T-Y*Z
42 E=A(1)*FND(B(2),C(2),B(3),C(3))
44 E=E-B(1)*FND(A(2),C(2),A(3),C(3))
46 E=E+C(1)*FND(A(2),B(2),A(3),B(3))
50 DISP “DET:”; E
52 IF E=0 THEN STOP
60 I(1,1)=FND(B(2),B(3),C(2),C(3))/E
62 I(1,2)=-FND(B(1),B(3),C(1),C(3))/E
64 I(1,3)=FND(B(1),B(2),C(1),C(2))/E
66 I(2,1)=-FND(A(2),A(3),C(2),C(3))/E
68 I(2,2)=FND(A(1),A(3),C(1),C(3))/E
70 I(2,3)=-FND(A(1),A(2),C(1),C(2))/E
72 I(3,1)=FND(A(2),A(3),B(2),B(3))/E
74 I(3,2)=-FND(A(1),A(3),B(1),B(3))/E
76 I(3,3)=FND(A(1),A(2),B(1),B(2))/E
78 IF H=2 THEN 100
80 FOR R=1 TO 3
82 FOR S=1 TO 3
84 DISP “I(“; R; “,”; S; “):”; I(R,S)
86 PAUSE
88 NEXT S
90 NEXT R
92 STOP
100 DIM Q(3)
102 FOR K=1 TO 3
104 Q(K)= I(K,1)*D(1) + I(K,2)*D(2) + I(K,3)*D(3)
106 DISP “Q”; K; “:”; Q(K) @ PAUSE
108 NEXT K

Example:

M = [[ 1, 2, -8 ] [ 0, -2, 9.5 ] [ 3.2, 2.7, -1 ]]
D = [[ 0.5 ] [ 1.5 ] [ 2.5 ]]

DET = -14.05
I ≈ [[ 1.6833, 1.3950, -0.2135 ] [ -2.1637, -1.7509, 0.6762 ] [ -0.4555, -0.2633, 0.1423 ]]

Solutions:
Q ≈ [[ 2.4004 ] [ -2.0178 ] [ -0.2669 ]]

Days From January 1, Days Between Dates

HP 71B – Day Counts

Number of Days from January 1

D = Day
M = Month
L = Leap Year indicator (1 if the year is a leap year, 0 if it is not)

Program DAYJAN1 (183 bytes)
10 DESTROY D,M,L,N
12 INPUT “MONTH:”; M
14 INPUT “DAY:”; D
16 INPUT “LEAP? (Y=1,N=0):”; L
20 IF M>2 THEN 28
21 REM M<=2
22 N=IP(30.6*(M+13))+D-429
24 GOTO 32
27 REM M>2
28 N=IP(30.6*(M+1))+D+L-64
32 DISP “# DAYS:”; N

Test 1: January 1 to May 29, non-leap year (L=0). Result: 148
Test 2: January 1 to May 29, leap year (L=1). Result: 149

Days between Dates

M, D, Y: Month, Day, four-digit year

Program DDAYS (299 bytes)
10 DESTROY M1, M2, D1, D2, Y1, Y2, F1, F2
11 DESTROY F,M,D,Y,N,X,Z
12 INPUT “1: M,D,Y:”; M1, D1, Y1
13 INPUT “2: M,D,Y:”; M2, D2, Y2
15 M=M1 @ D=D1 @ Y=Y1
17 GOSUB 40
19 F1=F
21 M=M2 @ D=D2 @ Y=Y2
23 GOSUB 40
25 F2=F
27 N=F2-F1
29 DISP “# DAYS:”; N
31 STOP
40 IF M>2 THEN X=IP(.4*M+2.3) ELSE X=0
42 IF M>2 THEN Z=Y ELSE Z=Y-1
44 F=365*Y+31*(M-1)+D+IP(Z/4)-X
46 RETURN

Test 1: January 2, 2015 to March 17, 2016 (1,2,2015 to 3,17,2016). Result: 440
Test 2: March 14, 1977 to June 29, 2016 (3,14,1977 to 6,29,2016). Result: 14352

Great Circle Distance

N: Longitude
E: Latitude

Separate Degrees (H), Minutes (M), and Seconds (S) during input

Program GRCIC (367 bytes)
10 DESTROY D,M,H,S
11 DESTROY N1,N2,E1,E2,G
12 DEGREES
14 DISP “N: LATITUDE” @ WAIT 1
16 DISP “E: LONGITUDE” @ WAIT 1
20 INPUT “N1: D,M,S:”; H,M,S
21 GOSUB 50
22 N1=D
25 INPUT “E1: D,M,S:”; H,M,S
26 GOSUB 50
27 E1=D
30 INPUT “N2: D,M,S:”; H,M,S
31 GOSUB 50
32 N2=D
35 INPUT “E2: D,M,S:”; H,M,S
36 GOSUB 50
37 E2=D
40 REM CALCULATION
41 G=SIN(N1)*SIN(N2)+COS(N1)*COS(N2)*COS(E1-E2)
42 G=ACOS(G)*3959*PI/180
44 DISP “DIST: “; G; “ MI”
46 STOP
50 D=SGN(H)*(ABS(H)+M/60+S/3600)
52 RETURN

Test:
Los Angeles, N = 34°13’0” and E = -118°15’0”
San Francisco, N = 37°47’0” and E = -112°25’0”
Distance ≈ 408.5961 mi

Euclid Division – Finding the GCD

Finds the GCD (greatest common divisor) between integers M and N.  The program displays a “calculating” screen while the calculation is in process.

Program EUCLID (143 bytes)

10 DESTROY M,N,A,B,C

15 INPUT “M,N:”; M,N

20 IF M>N THEN A=M  @ B=N

22 IF M<N THEN A=N @ B=M

30 C= A – IP(A/B)*B

35 DISP C;B;A   // this is the “busy” indicator

40 IF C=0 THEN 50

45 A=B @ B=C @ GOTO 30

50 DISP “GCD:”; B

Test 1:  M=144, N=14;  Result: 2

Test 2:  N=14, M=144;  Result: 2

Fractions:  Addition and Multiplication

Adds or multiplies two fractions W/X and Y/Z.  Gives the result in the simplest form.  Proper or improper fractions only.

Separate each part with a comma (W,X,Y,Z) as prompted.

Program FRAC (380 bytes)

10 DESTROY W,X,Y,Z,N,D,H,A,B,C

20 INPUT “1. + 2. *:”,H

24 ON H GOTO 41,51

41 REM ADD

42 INPUT “W/X+Y/Z:”; W,X,Y,Z

44 N=W*Z+X*Y @ D=X*Y

46 GOSUB 61

48 DISP “=”; N; “/”; D @ STOP

51 REM MULT

52 INPUT “W/X*Y/Z:”; W

54 N=W*Y @ D=X*Z

56 GOSUB 61

58 DISP “=”; N; “/”; D @ STOP

61 REM SIMPLIFY

62 IF N>D THEN A=N @ B=D

64 IF D>N THEN A=D @ B=N

66 IF D=N THEN N=1 @ D=1 @ RETURN

68 C= A – IP(A/B)*B

69 DISP C;B;A

70 IF C=0 THEN 74

72 A=B @ B=C @ GOTO 68

74 N=N/B @ D=D/B @ RETURN

Test 1:  4/7 + 3/13;   Input:  4,7,3,13.  Result:  73/91

Test 2: 4/7 * 3/13; Input: 4,7,3,13.  Result:  12/91

Statistics: Regression

The program CURVEFIT fits data to one of five regression models:

1.  Linear Regression:  y = a + bx

2.  Exponential Regression:  y = a * e^(b*x)

3.  Logarithm Regression: y = a + b * ln x

4.  Power Regression:  y = a * x^b

5.  Inverse Regression:  y = a + b/x

On the HP 71, the ln function is represented by LOG.

The program will allow different calculations with the same data set.

Program CURVEFIT (622 bytes)

10 DESTROY H,S,D,X,Y,A,B,E

12 STAT S(2) @ CLSTAT

20 REM CHOOSE REG

22 DISP “1. LIN 2. EXP 3. LOG” @ WAIT 1.5

24 DISP “4. POW 5. INV” @ WAIT 1.5

28 INPUT “CHOICE #:”; H

29 IF E=1 THEN 60

30 REM INPUT PROCESS

32 INPUT “X,Y:”; X,Y

34 ON H GOTO 36,38,40,42,44

36 ADD X,Y @ GOTO 46

38 ADD X,LOG(Y) @ GOTO 46

40 ADD LOG(X),Y @ GOTO 46

42 ADD LOG(X),LOG(Y) @ GOTO 46

44 ADD 1/X,Y @ GOTO 46

46 INPUT “DONE? (Y=1,N=0):”; D

48 IF D=0 THEN 32

60 LR 2,1,A,B

62 IF H=2 OR H=4 THEN A=EXP(A)

64 ON H GOTO 70,72,74,76,78

70 DISP A; “+”; B; “x” @ GOTO 80

72 DISP A; “*EXP(“; B; “x)” @ GOTO 80

74 DISP A; “+”; B; “*LOG(x)” @ GOTO 80

76 DISP A; “x^”; B @ GOTO 80

78 DISP A; “+”; B; “/x” @ GOTO 80

80 PAUSE

84 DISP “ANOTHER ANALYSIS?” @ WAIT 1.5

86 INPUT “Y=1,N=0:”; E

88 IF E=1 THEN 20

90 DISP “DONE”

 Have fun!  See you in the second half of 2016!  

Eddie

This blog is property of Edward Shore, 2016

Casio Classpad fx-CP400: Defining Functions

Casio Classpad fx-CP400:  Defining Functions

Using the Define command with the Classpad

A lot of formulas can be created using the Define command on the Casio Classpad.  For this blog I am using the fx-CP400.  My preference is to use the Main application, and call the Define command from the Action>Command submenu.  The expressions are stored as Functions.  If you want to transfer functions (and programs) from the fx-CP400 to a computer drive, export them to the Save-F folder first and then connect the Classpad to the computer.

Fibonacci Numbers

fibon(n) = approx((1+√5)^n – (1-√5)^n)/(2^n * √5)

The approx command is used in order to force a simplified numerical answer. 

Error Function

erf(z) = 2/√π * ∫ (e^-x^2 dx from 0 to z)

I think numerical integrals return approximate answers no matter what, please correct me if I am wrong.

Digital Root:  Counting all digits of a number, repeating the process until you get a single digit (0-9)

dr(n) = 1 + mod(n-1, 9)

I could only find the mod function in the catalog.  Unlike most Casio calculators, the Classpad’s mod function accepts negative numbers and non-integers.

Area of a Regular Polygon

aregpoly(n,s) = (n*s^2) / (4 * tan(180°/n))

The degree symbol (°) is needed to allow proper calculation regardless of angular mode.  The degree symbol is found in the Trig soft menu.

Great Circle: Distance between two places in kilometers

Note:

*  In order for the function to work properly, the Degree mode must be selected.

*  You can enter degrees in terms of degrees, minutes, seconds by choosing the DMS template.  This template is in the Math1 soft keyboard and is represented by three boxes (□ □ □) next to toDMS.  You can also call the template by pressing Action/Interactive, Transformation, DMS, dms.   The dms command’s syntax is  dms(degrees, minutes, seconds) and transforms the input into degree decimals.

Example:

Los Angeles,  N: 34°13’, E: -118°15’;   Tokyo,  N: 35°41’22.22”, E: 139°42’30.12”

Distance ≈ 8,803.002688 km

Formula (for km):

Distance ≈ cosˉ¹ (sin N1 sin N2 + cos N1 cos N2 cos (E1 – E2)) * 6371 * π / 180

If you want US miles, replace 6371 with 3959.

This blog is property of Edward Shore, 2016.

PrgCalcPro: Great Circle: Distance Between 2 Places on Earth in kilometers

PrgCalcPro: Great Circle: Distance Between 2 Places on Earth in kilometers

Important: Switch the calculator to Degrees mode before calculation! 

Formula:

E = acos( sin N1 sin N2 + cos N1 cos N2 cos (E1 - E2) )

Distance = 6371 * E * π / 180

Store the following values in decimal degrees:

A = N1 (longitude 1), B = E1 (latitude 1)

C = N2 (longitude 2), D = E2 (latitude 2)

Use [ K ] [ swap ] ( ° ' > to convert from H.M to Decimal degrees if needed)

Program:

   0: 6A  ;  RA

   1: 1C  ;  sin

   2: 6C  ;  RC

   3: 1C  ;  sin

   4: 12  ;  *

   5: 6A  ;  RA

   6: 1D  ;  cos

   7: 6C  ;  RC

   8: 1D  ;  cos

   9: 12  ;  *

  10: 6B  ;  RB

  11: 6D  ;  RD

  12: 11  ;  -

  13: 1D  ;  cos

  14: 12  ;  *

  15: 10  ;  +

  16: 1A  ;  acos

  17: 06  ;  6

  18: 03  ;  3

  19: 07  ;  7

  20: 01  ;  1

  21: 12  ;  *

  22: 20  ;  Pi

  23: 12  ;  *

  24: 01  ;  1

  25: 08  ;  8

  26: 00  ;  0

  27: 13  ;  /

  28: 50  ;  STOP

Example:

Los Angeles (N = 34.05°, E =-118.25°)

Moscow (N = 55.783333°, E = 37.616666°)

Approximate distance: 9,766.5757 km

This blog is property of Edward Shore, 2016. 

Casio fx-3650p: Programming - The Sequel

Casio fx-3650p: Programming - The Sequel

Time to revisit the Casio fx-3650p.  Since the last time I did this, Casio released an updated fx-3650p.  I still have the older version.  And still, the fx-3650p is not sold in the United States (sigh).   

Since the language of the fx-3650p is similar to the Casio graphing calculators and fx-5800p (in fact it is simplified, you can adopt and port these programs.)

Access to the first set of fx-3650P programs:

http://edspi31415.blogspot.com/2014/05/casio-fx-3650p-programming.html

Contents for this blog:

1.  Combination with Replacement

2.  Great Circle (Distance in km) 

3.  Orbital Speed and Period 

4.  Eccentricity and Area of an Ellipse

5.  Super Factorial

6.  Escape Velocity 

7.  Finance: Payment of a Monthly Mortgage

8.  Wind Chill Factor

9.  Speed of Sound in Dry Air 

Combination with Replacement 

Formula: nCr(A + B - 1, B)

Program: (17 steps)

? → A : ? → B :

(A + B - 1) [nCr] B 

Example: 

A = 17, B = 8, Result: 735,471

A = 52, B = 5, Result: 3,819,816

Great Circle (Distance in km)

Input:

A = latitude 1 (North is positive), B = longitude 1 (East is positive)

C = latitude 2, D = longitude 2

Program: (49 steps)

? → A : ? → B : ? → C : ? → D :

cos⁻¹ ( sin A sin C + cos A cos C cos (B - D → Y :

Y * 6371 * π ÷ 180 → Y

Example: Los Angels to Rome

Los Angeles:  Lat 34°03' N, Long 118°15' W (enter as negative)

Rome: Lat 41°54' N, Long 12°30' E

Result: (approx) 10,189.94397 km

Orbital Speed and Period   

Input:

A = Mass 1 (kg), B = Mass 2 (kg), D = Distance (m)

Output:

X = Orbital Speed (m/s), Y = Time for One Orbit to complete (in years)

Program: (61 steps)

? → A : ? → B : ? → D : 6.67384E-11 ( A + B → C :

√ ( C ÷ D → X ◢ 2 π √ ( D³ ÷ C → Y : Y ÷ 315576000 → Y

Example:

Sun: 1.989E30 kg (A)

Earth: 5.927E24 kg (B)

Avg. distance between Sun and Earth: 1.496E11 m (D)

Results:

Orbital Speed (Earth around Sun): 29787.91714 m/s (X)

Orbital Period: 0.999924841 yrs (Y)

Eccentricity and Area of an Ellipse

Ellipse with semi-axis lengths A and B, assuming that A ≤ B.  

Program: (25 steps)

? → A : ? → B : 

√ ( 1 - A² ÷ B²  → C ◢ π A B → D 

Example:

A = 5, B = 10

Results:

Eccentricity: 0.866025403 (C)

Area: 157.0796327 (D)

Super Factorial

Formula: spf(n) =   product(X!, X, 1, n)

For the Casio, n is an integer where 1 ≤ n ≤ 16

Source for formula:

Martin, Ángel M.  Sandmath_44: Math Extensions for the HP 41. Rev 44_E. 2012

Program: (35 steps)
? → A : 1 → Y : 1 → M : Lbl 0 :

M! * Y → Y : 1 M+ : A > M - 1 ⇒ Goto 0 : Y 

Examples:

spf(1) = 1

spf(2) = 2

spf(3) = 12

spf(4) = 288

Escape Velocity

Input:

A = mass (kg), D = radius (m)

Output:

Escape velocity in m/s

Program: (27 Steps)

? → A : ? → D : 

√ ( 2 * 6.67384E-11 * A ÷ D 

Example:

A = 5.927E24 kg (mass of Earth)

D = 6.371E3 m (radius of Earth)

Result: 11,143.37008 m/s

Finance: Payment of a Monthly Mortgage

Input:

A = Loan Amount

C = Annual Interest Rate 

D = Number of payments

Output:

B = Monthly Payment 

End of month payments are assumed.  No balloon amounts are assumed.  The program is set up to mimic behavior on most financial calculators (negative for outflows, positive for inflows).  

Program: (40 Steps)

? → A : ? → C : ? → D : C ÷ 1200 → X :

- A * X ÷ ( 1 - ( 1 + X ) ^ -D ) → B

Example:

A = 400,000

C = 5%

D = 360

Result:

B = -2,147.286491 (payment of 2,147.29 per month)

Wind Chill Factor (US Units)

Source for formula:

Glover, Thomas J.  Pocket Ref. 4th Edition.  Sequoia Publishing Inc.  Littleton, CO.  2012

Input:

A = temperature in °F

B = speed of the wind in mph

Output:  

C = Wind Chill in °F

Program: (47 steps)

? → A : ? → B : 35.74 + .6215 A 

- 35.75 B^.16 + .4275 * A * B^.16 → C 

Example:

A = 45°F

B = 15 mph

Result:

C = 38.23993448 °F

Speed of Sound in Dry Air

C = √ ( γ R (T + 273.15))

γ = 1.4 (ratio of specific heat of dry air)

R = 286.9 J/(kg K) (Individual Gas Constant)

T = temperature in °C

Sources:

"Speed of Sound" NASA.  Glenn Research Center. 

https://www.grc.nasa.gov/www/k-12/airplane/sound.html

Retrieved October 27, 2015

"The Individual and Universe Gas Constant".  The Engineering Toolbox. 

http://www.engineeringtoolbox.com/individual-universal-gas-constant-d_588.html

Retrieved October 27, 2015

Input:

C = temperature in °C

Output:

Speed of sound in m/s

Program: (25 Steps)

? → C : √(  1.4 * 286.9 * (C + 273.15

Example:

For 20°C, the speed of sound is approximately 343.1422868 m/s 

Enjoy!  Happy Halloween!  

Eddie 

This blog is property of Edward Shore. 2015