Home
JAQForum Ver 20.06
Log In or Join  
Active Topics
Local Time 01:25 21 Apr 2024 Privacy Policy
Jump to

Notice. New forum software under development. It's going to miss a few functions and look a bit ugly for a while, but I'm working on it full time now as the old forum was too unstable. Couple days, all good. If you notice any issues, please contact me.

Forum Index : Microcontroller and PC projects : DOS MMBasic Feature Request- Single Precision Math

Author Message
jwettroth

Regular Member

Joined: 02/08/2011
Location: United States
Posts: 70
Posted: 04:23pm 09 Apr 2021
Copy link to clipboard 
Print this post

I love DOS MMBasic and use it a lot- the ability to access serial ports is just great.  The ease of putting together little calculation programs quickly is awesome.

Degrading the accuracy of the floating point may seem like a strange request- I'm not trying to save memory, I am trying to develop a program on the PC and get an idea of how it will work when ported down to a Micromite.

I'm building a project that plays bird calls from 40 minutes before sunset to sunset.  For the prototype, I calculated the times off line by inputting my latitude and longitude from the web.  I then had an array of alarm times for each day during the play season- about 8 weeks-t worked fine but wasn't very elegant.

Over the winter, I did a little research and wrote a DOS MMBASIC app that will do all the trig calculations given your position.  It works well and agrees with the my web resources.  Its not very polished but I'd be happy to upload it.  I'm planning to port all of this back to a Micromite II for the final system. My DOS program and my Micromite program don't agree because the mite uses single precision. Would it be possible to have a data type or an option in DOS MMBASIC to use single precision floating point like the mites. Perhaps also, maybe there is a work around that I'm not seeing.  I was thinking about doing something like rounding the mantissa by scaling or similar.  Any assistance appreciated.

Regards All
John Wettroth
 
matherp
Guru

Joined: 11/12/2012
Location: United Kingdom
Posts: 8570
Posted: 05:15pm 09 Apr 2021
Copy link to clipboard 
Print this post

If you download the source from MMBASIC.COM and the free WATCOM C compiler you can rebuild the DOS version for single precision. It is just one change in the file configuration.h. Geoff has provided a batch file to do the build and it really is very easy.

#define MMFLOAT double                              // precision of all floating point operations
#define STR_SIG_DIGITS 9                            // number of significant digits to use when converting MMFLOAT to a string


becomes

#define MMFLOAT float                              // precision of all floating point operations
#define STR_SIG_DIGITS 5                            // number of significant digits to use when converting MMFLOAT to a string
 
jwettroth

Regular Member

Joined: 02/08/2011
Location: United States
Posts: 70
Posted: 06:08pm 09 Apr 2021
Copy link to clipboard 
Print this post

Thank you- that's a great solution.  I've always wanted to delve into compiling the source too.  Thanks

JW
John Wettroth
 
CaptainBoing

Guru

Joined: 07/09/2016
Location: United Kingdom
Posts: 1985
Posted: 07:18pm 09 Apr 2021
Copy link to clipboard 
Print this post

you might want to consider this: http://www.fruitoftheshed.com/MMBasic.Quick-and-Dirty-Daylight-Indicator.ashx

it uses key sun up/sundown times for your location  (stored in the case statements) then interpolates between the two based on the actual date.

It looks quite rough but is surprisingly accurate - to within five minutes or so - so much so you don't notice. I use it for numerous dusk-to-dawn light switches and I never have to bother with the lights - they just come on as twilight starts and go off as twilight ends in the morning.

Saves all the headache of a "proper" sunrise/set calculations using julian dates etc.

jus' sayin'
Edited 2021-04-10 05:19 by CaptainBoing
 
jwettroth

Regular Member

Joined: 02/08/2011
Location: United States
Posts: 70
Posted: 03:31am 10 Apr 2021
Copy link to clipboard 
Print this post

CaptainBoing-
Thanks I'll check it out.  The real calcuation is kind of ridiculous- got the algorithm from the US Naval Observatory.
John Wettroth
 
TassyJim

Guru

Joined: 07/08/2011
Location: Australia
Posts: 5886
Posted: 06:03am 10 Apr 2021
Copy link to clipboard 
Print this post

With most calculations there will be negligible difference between single an double precision.
Using this program, dating back to when MMBasic used line numbers.
5 DIM A(2), D(2)
10' Sunrise-Sunset adapted by TassyJim TBS
20 GOSUB 300
30 INPUT "Lat-deg="; B5
32 INPUT "Long-deg="; L5
40 INPUT "Time zone-hrs="; H
50 L5 = L5 / 360: Z0 = -H / 24
60 GOSUB 1170: T = (J - 2451545) + F
70 TT = T / 36525 + 1 ' TT = centuries
80 ' from 1900.0
90 GOSUB 410: T = T + Z0
100 '
110 ' Get Sun's Position
120 GOSUB 910: A(1) = A5: D(1) = D5
130 T = T + 1
140 GOSUB 910: A(2) = A5: D(2) = D5
150 IF A(2) < A(1) THEN A(2) = A(2) + P2
160 Z1 = DR * 90.833 ' Zenith dist.
170 S = SIN(B5 * DR): C = COS(B5 * DR)
180 Z = COS(Z1): M8 = 0: W8 = 0: PRINT
190 A0 = A(1): D0 = D(1)
200 DA = A(2) - A(1): DD = D(2) - D(1)
210 FOR C0 = 0 TO 23
220 P = (C0 + 1) / 24
230 A2 = A(1) + P * DA: D2 = D(1) + P * DD
240 GOSUB 490
250 A0 = A2: D0 = D2: V0 = V2
260 NEXT
270 GOSUB 820 ' Special msg?
280 END
290 '
300 ' Constants

320 P1 = 3.14159265
322 P2 = 2 * P1
330 DR = P1 / 180
332 K1 = 15 * DR * 1.0027379
340 Ss$ = "Sunset at "
350 sR$ = "Sunrise at "
360 M1t$ = "No sunrise this date"
370 M2t$ = "No sunset this date"
380 M3t$ = "Sun down all day"
390 M4t$ = "Sun up all day"
400 RETURN
410 ' LST at 0h zone time
420 T0 = T / 36525
430 S = 24110.5 + 8640184.8 * T0
440 S = S + 86636.6 * Z0 + 86400 * L5
450 S = S / 86400: S = S - INT(S)
460 T0 = S * 360 * DR
470 RETURN
480 '
490 ' Test an hour for an event
500 L0 = T0 + C0 * K1: L2 = L0 + K1
510 H0 = L0 - A0: H2 = L2 - A2
520 H1 = (H2 + H0) / 2 ' Hour angle,
530 D1 = (D2 + D0) / 2 ' declination,
540 ' at half hour
550 IF C0 > 0 THEN 570
560 V0 = S * SIN(D0) + C * COS(D0) * COS(H0) - Z
570 V2 = S * SIN(D2) + C * COS(D2) * COS(H2) - Z
580 IF SGN(V0) = SGN(V2) THEN 800
590 V1 = S * SIN(D1) + C * COS(D1) * COS(H1) - Z
600 Ax = 2 * V2 - 4 * V1 + 2 * V0: B = 4 * V1 - 3 * V0 - V2
610 Dy = B * B - 4 * Ax * V0: IF Dy < 0 THEN 800
620 Dy = SQR(Dy)
630 IF V0 < 0 AND V2 > 0 THEN PRINT sR$;
640 IF V0 < 0 AND V2 > 0 THEN M8 = 1
650 IF V0 > 0 AND V2 < 0 THEN PRINT Ss$;
660 IF V0 > 0 AND V2 < 0 THEN W8 = 1
670 E = (0-B + Dy) / (2 * Ax)
680 IF E > 1 OR E < 0 THEN E = (0-B - Dy) / (2 * Ax)
690 T3 = C0 + E + 1 / 120 ' Round off
700 H3 = INT(T3): M3 = INT((T3 - H3) * 60)
705  Hr$=STR$(H3):Hr$=RIGHT$("0"+Hr$,2):Mn$=STR$(M3):Mn$=RIGHT$(" 0"+Mn$,2)
710 PRINT Hr$;":";Mn$;
720 H7 = H0 + E * (H2 - H0)
730 N7 =0 -COS(D1) * SIN(H7)
740 D7 = C * SIN(D1) - S * COS(D1) * COS(H7)
750 AZ = ATN(N7 / D7) / DR
760 IF D7 < 0 THEN AZ = AZ + 180
770 IF AZ < 0 THEN AZ = AZ + 360
780 IF AZ > 360 THEN AZ = AZ - 360
782 PRINT ", azimuth ";
790 PRINT  AZ
800 RETURN
810 '
820 ' Special-message routine
830 IF M8 = 0 AND W8 = 0 THEN 870
840 IF M8 = 0 THEN PRINT M1t$
850 IF W8 = 0 THEN PRINT M2t$
860 GOTO 890
870 IF V2 < 0 THEN PRINT M3t$
880 IF V2 > 0 THEN PRINT M4t$
890 RETURN
900 '
910 ' Fundamental arguments
920 ' (Van Flandern &
930 ' Pulkkinen, 1979)
940 L = .779072 + .00273790931 * T
950 G = .993126 + .0027377785 * T
960 L = L - INT(L): G = G - INT(G)
970 L = L * P2: G = G * P2
980 V = .39785 * SIN(L)
990 V = V - .01 * SIN(L - G)
1000 V = V + .00333 * SIN(L + G)
1010 V = V - .00021 * TT * SIN(L)
1020 U = 1 - .03349 * COS(G)
1030 U = U - .00014 * COS(2 * L)
1040 U = U + .00008 * COS(L)
1050 W = -.0001 - .04129 * SIN(2 * L)
1060 W = W + .03211 * SIN(G)
1070 W = W + .00104 * SIN(2 * L - G)
1080 W = W - .00035 * SIN(2 * L + G)
1090 W = W - .00008 * TT * SIN(G)
1100 '
1110 ' Compute Sun's RA and Dec
1120 S = W / SQR(U - V * V)
1130 A5 = L + ATN(S / SQR(1 - S * S))
1140 S = V / SQR(U): D5 = ATN(S / SQR(1 - S * S))
1150 R5 = 1.00021 * SQR(U)
1160 RETURN
1165 '
1170 ' Calendar --> JD
1180 INPUT "Year= "; Y
1183 INPUT "Month= "; M
1186 INPUT "Day= "; Dy

1190 G = 1: IF Y < 1583 THEN G = 0
1200 D1 = INT(Dy): F = Dy - D1 - .5
1210 J =0 -INT(7 * (INT((M + 9) / 12) + Y) / 4)
1220 IF G = 0 THEN 1260
1230 S = SGN(M - 9): Ax = ABS(M - 9)
1240 J3 = INT(Y + S * INT(Ax / 7))
1250 J3 =0 -INT((INT(J3 / 100) + 1) * 3 / 4)
1260 J = J + INT(275 * M / 9) + D1 + G * J3
1270 J = J + 1721027 + 2 * G + 367 * Y
1280 IF F >= 0 THEN 1300
1290 F = F + 1: J = J - 1
1300 RETURN
1310 '
1320 ' This program by Roger W. Sinnott calculates the times of sunrise
1330 ' and sunset on any date, accurate to the minute within several
1340 ' centuries of the present. It correctly describes what happens in the
1350 ' arctic and antarctic regions, where the Sun may not rise or set on
1360 ' a given date. Enter north latitudes positive, west longitudes
1370 ' negative. For the time zone, enter the number of hours west of
1380 ' Greenwich (e.g., -5 for EST, -4 for EDT). The calculation is
1390 ' discussed in Sky & Telescope for August 1994, page 84.


The results are:
  Quote  Single
Sunrise at 06:41, azimuth  80.2146
Sunset at 17:54, azimuth  280.005

Double
Sunrise at 06:41, azimuth  80.21469347
Sunset at 17:54, azimuth  280.0046781


I really must get rid of the line numbers.

Jim
VK7JH
MMedit   MMBasic Help
 
matherp
Guru

Joined: 11/12/2012
Location: United Kingdom
Posts: 8570
Posted: 07:00am 10 Apr 2021
Copy link to clipboard 
Print this post

  Quote  With most calculations there will be negligible difference between single an double precision.


But not all  

Cdeagle's solar eclipse program gives no result at all with his reference example in single precision but works perfectly in double
 
jwettroth

Regular Member

Joined: 02/08/2011
Location: United States
Posts: 70
Posted: 04:16am 11 Apr 2021
Copy link to clipboard 
Print this post

Below is my code- heavily borrowed from sources in listing and pretty unpolished.

>>>>
'Sunrise/Sunset Algorithm
' john wettroth oct 2020
' algorithm and theory from Almanac for Computers USNO Washington DC 20392
' and heavily cribbed from edwilliams.org (great resource)


'  Written for MMBASIC for DOS-
'  consider early alpha- basic operation verified for my location but
'  no real testing done and no edge cases looked at.  Proceed at your own risk.
'
'Inputs
' day month year
'latitude longitude

outpu't
'zenith- type of sunrise/sunset
'official       90 deg 50'
'civil          96 deg
'nautical       102 deg
'astronomical   108 deg

'Notes
'  1. longitude is positive east
'  2. angles in degrees
'
'********************************************
'constants and inputs- these are hard coded here for testing
'could make them an input statement, etc
'
mon%=9 'october
day%=16 '16
year%=2020

'lat and log of location (Raleigh, NC)- positive is east, negative west
latitude=  35.7796
longitude= -78.6382

'PI/180 convert degrees to radians
PR=Pi/180


'BEGIN CODE
' day of year calculation
n1=Int(275* mon%/9)
n2=Int((mon%+9)/12)
n3=(1+Int((year%-4 * Int(year%/4) +2)/2))
n=n1-(n2*n3)+day% - 30
Print "Day of Year for ";Str$(mon%,2)+"/"+Str$(day%,2)+"/"+Right$(Str$(year%,4),2)+" is ";n

'convert our longitude to approximate time
lnghour=longitude /15
'rise time is N+ ((6-lnghour)/24)
'seTime
at= N+ ((18-lnghour)/24)
Print "Tvalue for Sunset =";at

'sun mean anomaly
M=(.9856*at)-3.289
Print "Mean Solar Anomaly=";M


'sun true longitude
L= M+ (1.916 * Sin(M*PR)) + (.020 * Sin(2*M*PR)) + 282.634
L=L-(360*(L\360))

Print "sun true longitude=";L

'sun right ascension
RA=(180/Pi)*Atn(.91764 * Tan(pr*L))
Print "right ascension=";ra

'ra adjusted to l quadrant
LQUAD= 90* Int(L/90): Print "LQUAD =";LQUAD
RQUAD= 90* Int(RA/90):Print "RQUAD =";RQUAD
RA=RA+(LQUAD-RQUAD)
Print "RA Adjusted=";RA

'CONVERT RA TO HOURS
RA=RA/15
Print "RA IN HOURS";RA

'DECLINATION
SDEC= .39782*Sin(L*PR)
CDEC= Cos(pr*(ASin(SDEC)/pr))
Print "SINDEC =";SDEC
Print "COSDEC =";CDEC

'SUN LOCAL HOUR ANGLE
COSH = (Cos(90.833*PR) - (SDEC*Sin(LATITUDE*PR)))/(CDEC*Cos(LATITUDE*PR))
If COSH > 1 Then Print "SUN NEVER RISES AT THIS LOCATION"
If COSH < -1 Then Print "SUN NEVER SETS HERE"
Print "Cos of Local Hour Angle =";cosh

'convert cosh to hours
H=180*ACos(cosh)/Pi
Print "Hour Angle in Degrees =";H
H=H/15
Print "Hour angle in time =";H
'local mean time
T=H+RA-(.06571*at)-6.622
If t>24 Then t=t-24
If t<0 Then T=t+24
Print "SUNSET LOCAL MEAN TIME=";T

'adj to utc and subtract lnghour
T=T-LNGHOUR
Print "UTC TIME FOR SUNSET";T
LT=T-5+1
Print "TIME OF SUNSET IS ";Int(LT);" HOURS AND "; Int(60*(LT-Int(LT))); " MINUTES"

End
John Wettroth
 
Print this page


To reply to this topic, you need to log in.

© JAQ Software 2024