Error 47 undefined label

[DaveJJ]

Hi,

After many years away from programming I decided to return to it using Freebasic as I sort of hoped that it would compile QB45 code relatively easily.
I chose a program which predicts sunrise/sunset time in QuickBasic for a project I fancy on the Raspberry

After much editing I came to the undefined label error which I don't understand and can't find explained simply elsewhere. As i understand it a label is the name of a function or sub etc but the erro refers to a variable as follows

FUNCTION cn(x AS DOUBLE) as double
dim z as double
z = Cos(x * 0.0174532925199433)
return z
END FUNCTION

giving compiler output

/home/pi/Freebasic/riset.bas(177) warning 13(0): Function result was not explicitly set
/home/pi/Freebasic/riset.bas() error 47: Undefined label, Z

Which don't make sense to me, can anyone help explain what the problem really is.

Thanks

[dkl]

Hello,

It's certainly a weird error message, but perhaps I can clear it up.

In FreeBASIC's -lang qb mode, Return is not supported for returning function results. Instead, Return is the counter-part to Gosub, and function results have to be assigned as in QB:

Code: Select all

myFunction = resultvalue

This seems to work:

Code: Select all

#lang "qb"
FUNCTION cn(x AS DOUBLE) as double
	dim z as double
	z = Cos(x * 0.0174532925199433)
	cn = z
END FUNCTION

[DaveJJ]

Superquick reply.

Many thanks. I started to come to the same conclusion and altered the script BACK to funcname=xxx from the Return xxxx that it was and the errors are altering. I'll finish the task and report back later ... food calls. Shame I tidied the original up so much as I came across the -lang qb switch a little too late to rescue me from my enthusiasm

I will return.

[DaveJJ]

Well, I did all the editing I could see and reduced the errors to these few (no doubt they may increase if I solve the problem with these)

the compiler errors are now no longer too many to deal with but remain as

Code: Select all

/home/pi/Freebasic/riset.bas(13) error 4: Duplicated definition
DECLARE SUB quad (ym AS DOUBLE, y0 AS DOUBLE, yp AS DOUBLE, xe AS DOUBLE, ye AS DOUBLE, z1 AS DOUBLE, z2 AS DOUBLE, nz AS INTEGER)
                                                                                                                                 ^
/home/pi/Freebasic/riset.bas(212) error 4: Duplicated definition, at parameter 1 (mjd) of LMST()
FUNCTION lmst(mjd AS DOUBLE, glong AS DOUBLE) as double
                                                 ^
/home/pi/Freebasic/riset.bas(217) error 1: Argument count mismatch, found ')'
mjd0 = ipart(mjd)
                ^
/home/pi/Freebasic/riset.bas(218) error 1: Argument count mismatch, found '-'
ut = (mjd - mjd0) * 24
          ^

the code is very long to look at but I include it because the inclusion of symbols ! $ # % & after variables and constants confuses me. Here is the text

Code: Select all

DECLARE FUNCTION hm(ut AS DOUBLE) as double
DECLARE FUNCTION sinalt(iobj AS INTEGER, mjd0 AS DOUBLE, hour AS DOUBLE, glong AS DOUBLE, cphi AS DOUBLE, sphi AS DOUBLE) as double
DECLARE SUB quad (ym AS DOUBLE, y0 AS DOUBLE, yp AS DOUBLE, xe AS DOUBLE, ye AS DOUBLE, z1 AS DOUBLE, z2 AS DOUBLE, nz AS INTEGER)
DECLARE FUNCTION fpart (x AS DOUBLE) as double
DECLARE FUNCTION lmst(mjd AS DOUBLE, lambda AS DOUBLE) as double
DECLARE FUNCTION calday(mjd AS DOUBLE) as string
DECLARE FUNCTION ipart(x AS DOUBLE) as double
DECLARE FUNCTION cn(x AS DOUBLE) as double
DECLARE FUNCTION mjd(y AS INTEGER, m AS INTEGER, d AS INTEGER, h AS DOUBLE) as double
DECLARE FUNCTION sn(x AS DOUBLE) as double
DECLARE SUB moon (t AS DOUBLE, ra AS DOUBLE, dec AS DOUBLE)
DECLARE SUB sun (t AS DOUBLE, ra AS DOUBLE, dec AS DOUBLE)
DECLARE SUB quad (ym AS DOUBLE, y0 AS DOUBLE, yp AS DOUBLE, xe AS DOUBLE, ye AS DOUBLE, z1 AS DOUBLE, z2 AS DOUBLE, nz AS INTEGER)
'       Rise and set times for Sun and Moon
'       Adapted and modified from Montenbruck
'       and Pfleger, 'Astronomy on the personal
'       Computer' 3rd Edition, Springer
'       section 3.8
'       Accuracy of detection of 'always below' and 'always above'
'       situations depends on the approximate routines used for Sun
'       and Moon. For instance, 1999 Dec 25th, at 0 long, 67.43 lat
'       this program will give an 8 minute long day between sunrise
'       and sunset. More accurate programs say the Sun is always below
'       the horizon on this day.
'
dim as integer y,m,d

p$ = "    ####"
DEFDBL A-Z
pi = 4 * ATN(1)
rads = pi / 180
degs = 180 / pi
DIM sinho(3)
DIM obname$(5)
obname$(1) = "Moon"
obname$(2) = "Sun"
obname$(3) = "Nautical twilight"
CLS
PRINT "   Rise and set for Sun and Moon"
PRINT "   ============================="
PRINT
INPUT "   Year (yyyy) - - - - - - - - :", y%
INPUT "   Month  (mm) - - - - - - - - :", m%
INPUT "   Day    (dd) - - - - - - - - :", d%
INPUT "   Time zone (East +) - - - -  :", zone
INPUT "   Longitude (w neg, decimals) :", glong
INPUT "   Latitude  (n pos, decimals) :", glat
glong = -glong  'routines use east longitude negative convention
zone = zone / 24
date = mjd(y, m, d, 0) - zone
'define the altitudes for each object
'treat twilight as a separate object 3, so sinalt routine
'falls through to finding Sun altitude again
sl = sn(glat)
cl = cn(glat)
sinho(1) = sn(8! / 60!)         'moonrise - average diameter used
sinho(2) = sn(-50! / 60!)       'sunrise - classic value for refraction
sinho(3) = sn(-12!)             'nautical twilight
xe = 0
ye = 0
z1 = 0
z2 = 0
FOR iobj% = 1 TO 3
    utrise = 0
    utset = 0
    rise = 0
    sett = 0
    hour = 1
    zero2 = 0
    ' See STEP 1 and 2 of Web page description.
    ym = sinalt(iobj%, date, hour - 1, glong, cl, sl) - sinho(iobj%)
    IF ym > 0! THEN above = 1 ELSE above = 0
    'used later to classify non-risings
    DO
        'STEP 1 and STEP 3 of Web page description 
        y0 = sinalt(iobj%, date, hour, glong, cl, sl) - sinho(iobj%)
        yp = sinalt(iobj%, date, hour + 1, glong, cl, sl) - sinho(iobj%)
        xe = 0
        ye = 0
        z1 = 0
        z2 = 0
        nz% = 0
        'STEP 4 of web page description
        quad (ym, y0, yp, xe, ye, z1, z2, nz%)
        SELECT CASE nz%
            'cases depend on values of discriminant - inner part of STEP 4
            CASE 0 'nothing  - go to next time slot
            CASE 1                      ' simple rise / set event
                IF (ym < 0!) THEN       ' must be a rising event
                        utrise = hour + z1
                        rise = 1
                ELSE                    ' must be setting
                        utset = hour + z1
                        sett = 1
                END IF
            CASE 2                      ' rises and sets within interval
                IF (ye < 0!) THEN       ' minimum - so set then rise
                        utrise = hour + z2
                        utset = hour + z1
                ELSE                    ' maximum - so rise then set
                        utrise = hour + z1
                        utset = hour + z2
                END IF
                rise = 1
                sett = 1
                zero2 = 1
            END SELECT
        ym = yp     'reuse the ordinate in the next interval
        hour = hour + 2
    ' STEP 5 of Web page description - have we finished for this object?
    LOOP UNTIL (hour = 25) OR (rise * sett = 1)
    utrise = hm(utrise)
    utset = hm(utset)
    'STEP 6 of Web page description
    PRINT
    PRINT "   "; obname$(iobj%)
    ' logic to sort the various rise and set states
    IF (rise = 1 OR sett = 1) THEN   'current object rises and sets today
        IF rise = 1 THEN
            PRINT USING p$; utrise
        ELSE
            PRINT "    ----"
        END IF
        IF sett = 1 THEN
            PRINT USING p$; utset
        ELSE
            PRINT "    ----"
        END IF
    ELSE              'current object not so simple
        IF above = 1 THEN
            SELECT CASE iobj%
                CASE 1, 2: PRINT "    always above horizon"
                CASE 3: PRINT "    always bright"
            END SELECT
        ELSE
            SELECT CASE iobj%
                CASE 1, 2: PRINT "    always below horizon"
                CASE 3: PRINT "    always dark"
            END SELECT
        END IF
    END IF
'STEP 7 of Web page description
NEXT iobj%
END
                             

DEFSNG A-Z
FUNCTION calday(x AS DOUBLE) as string
'    returns calendar date as a string in international format
'    given the modified julian date
'    BC dates are in calendar format - i.e. no year zero
'    Gregorian dates are returned after 1582 Oct 10th
'    In English colonies and Sweeden, this does not reflect
'    historical dates
dim dt as string
jd# = x + 2400000.5#
jd0 = ipart(jd# + .5)
IF jd0 < 2299161# THEN
    c = jd0 + 1524#
ELSE
    b = ipart((jd0 - 1867216.25#) / 36524.25#)
    c = jd0 + (b - ipart(b / 4)) + 1525#
END IF
d = ipart((c - 122.1#) / 365.25#)
e = 365# * d + ipart(d / 4)
F = ipart((c - e) / 30.6001)
day = ipart(c - e + .5) - ipart(30.6001 * F)
month = F - 1 - 12 * ipart(F / 14)
year = d - 4715 - ipart((month + 7) / 10)
dt = STR$(year) + STR$(month) + STR$(day)
calday = dt
'return dt
END FUNCTION

FUNCTION cn(x AS DOUBLE) as double
dim z as double

z = Cos(x * 0.0174532925199433)
cn = z
'return z
END FUNCTION

DEFDBL A-Z
FUNCTION fpart (x AS DOUBLE) as double
'       returns fractional part of a number
x = x - INT(x)
IF x < 0 THEN
   x = x + 1
END IF
fpart= x
END FUNCTION

FUNCTION hm (ut AS DOUBLE) as double
' returns number containing the time written in hours and minutes
' rounded to the nearest minute
dim as double h, m

ut = Int(ut * 60! + .5) / 60!   'round ut to nearest minute
h = Int(ut)
m = Int(60! * (ut - h) + .5)
hm=int(100*h+m)
'return Int(100 * h + m)
END FUNCTION

DEFSNG A-Z
FUNCTION ipart(x AS DOUBLE) as double
ipart=sgn(x) * Int(Abs(x))
'return SGN(x) * INT(ABS(x))
END FUNCTION

DEFDBL A-Z
FUNCTION lmst(mjd AS DOUBLE, glong AS DOUBLE) as double
'    returns the local siderial time for
'    the mjd and longitude specified
dim as double mjd0

mjd0 = ipart(mjd)
ut = (mjd - mjd0) * 24
t = (mjd0 - 51544.5) / 36525
gmst = 6.697374558 + 1.0027379093# * ut
gmst = gmst + (8640184.812866 + (0.093104 - 0.0000062 * t) * t) * t / 3600#
lmst= (24 * fpart((gmst - glong / 15) / 240))
'return 24 * fpart((gmst - glong / 15) / 24)
END FUNCTION

DEFSNG A-Z
FUNCTION mjd(y AS INTEGER, m AS INTEGER, d AS INTEGER, h AS DOUBLE) as double
'   returns modified julian date
'   number of days since 1858 Nov 17 00:00h
'   valid for any date since 4713 BC
'   assumes gregorian calendar after 1582 Oct 15, Julian before
'   Years BC assumed in calendar format, i.e. the year before 1 AD is 1 BC
a# = 10000# * y + 100# * m + d
IF y < 0 THEN y = y + 1
IF m <= 2 THEN
   m = m + 12
   y = y - 1
END IF
IF a# <= 15821004.1# THEN
   b = -2 + (y + 4716) \ 4 - 1179
ELSE
   b = (y \ 400) - (y \ 100) + (y \ 4)
END IF
a = 365 * y - 679004
mjd = (a + b + ipart(30.6001 * (m + 1)) + d + h / 24)
'return a# + b + ipart(30.6001# * (m + 1)) + d + h / 24
END FUNCTION

DEFDBL A-Z
SUB moon (t AS DOUBLE, ra AS DOUBLE, dec AS DOUBLE)
' returns ra and dec of Moon to 5 arc min (ra) and 1 arc min (dec)
' for a few centuries either side of J2000.0
' Predicts rise and set times to within minutes for about 500 years
' in past - TDT and UT time diference may become significant for long
' times
p2 = 6.283185307#
ARC = 206264.8062#
COSEPS = .91748
SINEPS = .39778
L0 = fpart(.606433 + 1336.855225# * t)    'mean long Moon in revs
L = p2 * fpart(.374897 + 1325.55241# * t) 'mean anomaly of Moon
LS = p2 * fpart(.993133 + 99.997361# * t) 'mean anomaly of Sun
d = p2 * fpart(.827361 + 1236.853086# * t)'diff longitude sun and moon
F = p2 * fpart(.259086 + 1342.227825# * t)'mean arg latitude
' longitude correction terms
dL = 22640 * SIN(L) - 4586 * SIN(L - 2 * d)
dL = dL + 2370 * SIN(2 * d) + 769 * SIN(2 * L)
dL = dL - 668 * SIN(LS) - 412 * SIN(2 * F)
dL = dL - 212 * SIN(2 * L - 2 * d) - 206 * SIN(L + LS - 2 * d)
dL = dL + 192 * SIN(L + 2 * d) - 165 * SIN(LS - 2 * d)
dL = dL - 125 * SIN(d) - 110 * SIN(L + LS)
dL = dL + 148 * SIN(L - LS) - 55 * SIN(2 * F - 2 * d)
' latitude arguments
S = F + (dL + 412 * SIN(2 * F) + 541 * SIN(LS)) / ARC
h = F - 2 * d
' latitude correction terms
N = -526 * SIN(h) + 44 * SIN(L + h) - 31 * SIN(h - L) - 23 * SIN(LS + h)
N = N + 11 * SIN(h - LS) - 25 * SIN(F - 2 * L) + 21 * SIN(F - L)
lmoon = p2 * fpart(L0 + dL / 1296000#)  'Lat in rads
bmoon = (18520# * SIN(S) + N) / ARC     'long in rads
' convert to equatorial coords using a fixed ecliptic
CB = COS(bmoon)
x = CB * COS(lmoon)
V = CB * SIN(lmoon)
W = SIN(bmoon)
y = COSEPS * V - SINEPS * W
Z = SINEPS * V + COSEPS * W
rho = SQR(1# - Z * Z)
dec = (360# / p2) * ATN(Z / rho)
ra = (48# / p2) * ATN(y / (x + rho))
IF ra < 0 THEN
        ra = ra + 24#
END IF
END SUB

SUB quad (ym AS DOUBLE, y0 AS DOUBLE, yp AS DOUBLE, xe AS DOUBLE, ye AS DOUBLE, z1 AS DOUBLE, z2 AS DOUBLE, nz AS INTEGER)
'  finds a parabola through three points and returns values of
'  coordinates of extreme value (xe, ye) and zeros if any (z1, z2)
'  assumes that the x values are -1, 0, +1
nz = 0
a = .5 * (ym + yp) - y0
b = .5 * (yp - ym)
c = y0
xe = -b / (2! * a)              'x coord of symmetry line
ye = (a * xe + b) * xe + c      'extreme value for y in interval
dis = b * b - 4! * a * c        'discriminant
IF dis > 0 THEN                 'there are zeros
    dx = .5 * SQR(dis) / ABS(a)
    z1 = xe - dx
    z2 = xe + dx
    IF (ABS(z1) <= 1!) THEN nz = nz + 1     'This zero is in interval
    IF (ABS(z2) <= 1!) THEN nz = nz + 1     'This zero is in interval
    IF (z1 < -1!) THEN z1 = z2
END IF
END SUB

FUNCTION sinalt (iobj AS INTEGER, mjd0 AS DOUBLE, hour AS DOUBLE, glong AS DOUBLE, cphi AS DOUBLE, sphi AS DOUBLE) as double
' returns sine of the altitude of either the sun or the moon given the
' modified julian day number at midnight UT and the hour of the UT day,
' the longitude of the observer, and the sine and cosine of the latitude
' of the observer
ra = 0
dec = 0
instant = mjd0 + hour / 24#
t = (instant - 51544.5#) / 36525#
IF (iobj = 1) THEN
        moon (t, ra, dec)
ELSE
        sun (t, ra, dec)
END IF
tau = 15# * (lmst(instant, glong) - ra)   'hour angle of object
sinalt= sphi * sn(dec) + cphi * cn(dec) * cn(tau)
'return sphi * sn(dec) + cphi * cn(dec) * cn(tau)
END FUNCTION

DEFSNG A-Z
FUNCTION sn(x AS DOUBLE) as double
sn=Sin(x*0.0174532925199433)
'return SIN(x * .0174532925199433#)
END FUNCTION

DEFDBL A-Z
SUB sun (t AS DOUBLE, ra AS DOUBLE, dec AS DOUBLE)
' Returns RA and DEC of Sun to roughly 1 arcmin for few hundred
' years either side of J2000.0
p2 = 6.283185307#
COSEPS = .91748
SINEPS = .39778
m = p2 * fpart(.993133 + 99.997361# * t)        'Mean anomaly
dL = 6893# * SIN(m) + 72# * SIN(2 * m)          'Eq centre
L = p2 * fpart(.7859453# + m / p2 + (6191.2# * t + dL) / 1296000#)
' convert to RA and DEC - ecliptic latitude of Sun taken as zero
sl = SIN(L)
x = COS(L)
y = COSEPS * sl
Z = SINEPS * sl
rho = SQR(1# - Z * Z)
dec = (360# / p2) * ATN(Z / rho)
ra = (48# / p2) * ATN(y / (x + rho))
IF ra < 0 THEN ra = ra + 24
END SUB

There are many constants which commence .xyzetc which I have started to change to 0.xyzetc as there seem to be compiler errors associated with them. Is this correct?

Many thanks

[DaveJJ]

PS I forgot to mention I don't understand the repeated

DEFSNG A-Z
and
DEFDBL A-Z

which occurs in various places, are they global definitions?

[dodicat]

DEFSNG A-Z

All variables starting with A to Z are singles.

I've got your code running, a couple of errors, you had quad declared twice, and you use mjd for different things.
Here's the runner (In code block form)

Code: Select all

#lang "qb"
DECLARE FUNCTION hm# (ut AS DOUBLE)
DECLARE FUNCTION sinalt# (iobj AS INTEGER, mjd0 AS DOUBLE, hour AS DOUBLE, glong AS DOUBLE, cphi AS DOUBLE, sphi AS DOUBLE)
'DECLARE SUB quad (ym AS DOUBLE, y0 AS DOUBLE, yp AS DOUBLE, xe AS DOUBLE, ye AS DOUBLE, z1 AS DOUBLE, z2 AS DOUBLE, nz AS INTEGER)
DECLARE FUNCTION fpart# (x AS DOUBLE)
DECLARE FUNCTION lmst# (_mjd AS DOUBLE, lambda AS DOUBLE)
DECLARE FUNCTION calday$ (mjd AS DOUBLE)
DECLARE FUNCTION ipart# (x AS DOUBLE)
DECLARE FUNCTION cn# (x AS DOUBLE)
DECLARE FUNCTION mjd# (y AS INTEGER, m AS INTEGER, d AS INTEGER, h AS DOUBLE)
DECLARE FUNCTION sn# (x AS DOUBLE)
DECLARE SUB moon (t AS DOUBLE, ra AS DOUBLE, dec AS DOUBLE)
DECLARE SUB sun (t AS DOUBLE, ra AS DOUBLE, dec AS DOUBLE)
DECLARE SUB quad (ym AS DOUBLE, y0 AS DOUBLE, yp AS DOUBLE, xe AS DOUBLE, ye AS DOUBLE, z1 AS DOUBLE, z2 AS DOUBLE, nz AS INTEGER)
' Rise and set times for Sun and Moon
' Adapted and modified from Montenbruck
' and Pfleger, 'Astronomy on the personal
' Computer' 3rd Edition, Springer
' section 3.8
' Accuracy of detection of 'always below' and 'always above'
' situations depends on the approximate routines used for Sun
' and Moon. For instance, 1999 Dec 25th, at 0 long, 67.43 lat
' this program will give an 8 minute long day between sunrise
' and sunset. More accurate programs say the Sun is always below
' the horizon on this day.
'
p$ = " ####"
DEFDBL A-Z
pi = 4 * ATN(1)
rads = pi / 180
degs = 180 / pi
DIM sinho(3)
DIM obname$(5)
obname$(1) = "Moon"
obname$(2) = "Sun"
obname$(3) = "Nautical twilight"
CLS
PRINT " Rise and set for Sun and Moon"
PRINT " ============================="
PRINT
INPUT " Year (yyyy) - - - - - - - - :", y%
INPUT " Month (mm) - - - - - - - - :", m%
INPUT " Day (dd) - - - - - - - - :", d%
INPUT " Time zone (East +) - - - - :", zone
INPUT " Longitude (w neg, decimals) :", glong
INPUT " Latitude (n pos, decimals) :", glat
glong = -glong 'routines use east longitude negative convention
zone = zone / 24
date = mjd(y%, m%, d%, 0#) - zone
'define the altitudes for each object
'treat twilight as a separate object 3, so sinalt routine
'falls through to finding Sun altitude again
sl = sn(glat)
cl = cn(glat)
sinho(1) = sn(8! / 60!) 'moonrise - average diameter used
sinho(2) = sn(-50! / 60!) 'sunrise - classic value for refraction
sinho(3) = sn(-12!) 'nautical twilight
xe = 0
ye = 0
z1 = 0
z2 = 0
FOR iobj% = 1 TO 3
utrise = 0
utset = 0
rise = 0
sett = 0
hour = 1
zero2 = 0
' See STEP 1 and 2 of Web page description.
ym = sinalt(iobj%, date, hour - 1, glong, cl, sl) - sinho(iobj%)
IF ym > 0! THEN above = 1 ELSE above = 0
'used later to classify non-risings
DO
'STEP 1 and STEP 3 of Web page description 
y0 = sinalt(iobj%, date, hour, glong, cl, sl) - sinho(iobj%)
yp = sinalt(iobj%, date, hour + 1, glong, cl, sl) - sinho(iobj%)
xe = 0
ye = 0
z1 = 0
z2 = 0
nz% = 0
'STEP 4 of web page description
quad ym, y0, yp, xe, ye, z1, z2, nz%
SELECT CASE nz%
'cases depend on values of discriminant - inner part of STEP 4
CASE 0 'nothing - go to next time slot
CASE 1 ' simple rise / set event
IF (ym < 0!) THEN ' must be a rising event
utrise = hour + z1
rise = 1
ELSE ' must be setting
utset = hour + z1
sett = 1
END IF
CASE 2 ' rises and sets within interval
IF (ye < 0!) THEN ' minimum - so set then rise
utrise = hour + z2
utset = hour + z1
ELSE ' maximum - so rise then set
utrise = hour + z1
utset = hour + z2
END IF
rise = 1
sett = 1
zero2 = 1
END SELECT
ym = yp 'reuse the ordinate in the next interval
hour = hour + 2
' STEP 5 of Web page description - have we finished for this object?
LOOP UNTIL (hour = 25) OR (rise * sett = 1)
utrise = hm(utrise)
utset = hm(utset)
'STEP 6 of Web page description
PRINT
PRINT " "; obname$(iobj%)
' logic to sort the various rise and set states
IF (rise = 1 OR sett = 1) THEN 'current object rises and sets today
IF rise = 1 THEN
PRINT USING p$; utrise
ELSE
PRINT " ----"
END IF
IF sett = 1 THEN
PRINT USING p$; utset
ELSE
PRINT " ----"
END IF
ELSE 'current object not so simple
IF above = 1 THEN
SELECT CASE iobj%
CASE 1, 2: PRINT " always above horizon"
CASE 3: PRINT " always bright"
END SELECT
ELSE
SELECT CASE iobj%
CASE 1, 2: PRINT " always below horizon"
CASE 3: PRINT " always dark"
END SELECT
END IF
END IF
'STEP 7 of Web page description
NEXT iobj%
sleep
END


DEFSNG A-Z
FUNCTION calday$ (x AS DOUBLE)
' returns calendar date as a string in international format
' given the modified julian date
' BC dates are in calendar format - i.e. no year zero
' Gregorian dates are returned after 1582 Oct 10th
' In English colonies and Sweeden, this does not reflect
' historical dates
jd# = x + 2400000.5#
jd0 = ipart(jd# + .5)
IF jd0 < 2299161# THEN
c = jd0 + 1524#
ELSE
b = ipart((jd0 - 1867216.25#) / 36524.25#)
c = jd0 + (b - ipart(b / 4)) + 1525#
END IF
d = ipart((c - 122.1#) / 365.25#)
e = 365# * d + ipart(d / 4)
F = ipart((c - e) / 30.6001)
day = ipart(c - e + .5) - ipart(30.6001 * F)
month = F - 1 - 12 * ipart(F / 14)
year = d - 4715 - ipart((month + 7) / 10)
calday$ = STR$(year) + STR$(month) + STR$(day)
END FUNCTION

FUNCTION cn# (x AS DOUBLE)
cn = COS(x * .0174532925199433#)
END FUNCTION

DEFDBL A-Z
FUNCTION fpart# (x AS DOUBLE)
' returns fractional part of a number
x = x - INT(x)
IF x < 0 THEN
x = x + 1
END IF
fpart = x
END FUNCTION

FUNCTION hm (ut AS DOUBLE)
' returns number containing the time written in hours and minutes
' rounded to the nearest minute
ut = INT(ut * 60! + .5) / 60! 'round ut to nearest minute
h = INT(ut)
m = INT(60! * (ut - h) + .5)
hm = INT(100 * h + m)
END FUNCTION

DEFSNG A-Z
FUNCTION ipart# (x AS DOUBLE)
ipart = SGN(x) * INT(ABS(x))
END FUNCTION

DEFDBL A-Z
FUNCTION lmst# (_mjd AS DOUBLE, glong AS DOUBLE)
' returns the local siderial time for
' the mjd and longitude specified
mjd0 = ipart(_mjd)
ut = (_mjd - mjd0) * 24
t = (mjd0 - 51544.5) / 36525
gmst = 6.697374558# + 1.0027379093# * ut
gmst = gmst + (8640184.812866# + (.093104 - .0000062 * t) * t) * t / 3600#
lmst = 24# * fpart((gmst - glong / 15#) / 24#)
END FUNCTION

DEFSNG A-Z
FUNCTION mjd# (y AS INTEGER, m AS INTEGER, d AS INTEGER, h AS DOUBLE)
' returns modified julian date
' number of days since 1858 Nov 17 00:00h
' valid for any date since 4713 BC
' assumes gregorian calendar after 1582 Oct 15, Julian before
' Years BC assumed in calendar format, i.e. the year before 1 AD is 1 BC
a# = 10000# * y + 100# * m + d
IF y < 0 THEN y = y + 1
IF m <= 2 THEN
m = m + 12
y = y - 1
END IF
IF a# <= 15821004.1# THEN
b = -2 + (y + 4716) \ 4 - 1179
ELSE
b = (y \ 400) - (y \ 100) + (y \ 4)
END IF
a# = 365# * y - 679004#
mjd = a# + b + ipart(30.6001# * (m + 1)) + d + h / 24
END FUNCTION

DEFDBL A-Z
SUB moon (t AS DOUBLE, ra AS DOUBLE, dec AS DOUBLE)
' returns ra and dec of Moon to 5 arc min (ra) and 1 arc min (dec)
' for a few centuries either side of J2000.0
' Predicts rise and set times to within minutes for about 500 years
' in past - TDT and UT time diference may become significant for long
' times
p2 = 6.283185307#
ARC = 206264.8062#
COSEPS = .91748
SINEPS = .39778
L0 = fpart(.606433 + 1336.855225# * t) 'mean long Moon in revs
L = p2 * fpart(.374897 + 1325.55241# * t) 'mean anomaly of Moon
LS = p2 * fpart(.993133 + 99.997361# * t) 'mean anomaly of Sun
d = p2 * fpart(.827361 + 1236.853086# * t)'diff longitude sun and moon
F = p2 * fpart(.259086 + 1342.227825# * t)'mean arg latitude
' longitude correction terms
dL = 22640 * SIN(L) - 4586 * SIN(L - 2 * d)
dL = dL + 2370 * SIN(2 * d) + 769 * SIN(2 * L)
dL = dL - 668 * SIN(LS) - 412 * SIN(2 * F)
dL = dL - 212 * SIN(2 * L - 2 * d) - 206 * SIN(L + LS - 2 * d)
dL = dL + 192 * SIN(L + 2 * d) - 165 * SIN(LS - 2 * d)
dL = dL - 125 * SIN(d) - 110 * SIN(L + LS)
dL = dL + 148 * SIN(L - LS) - 55 * SIN(2 * F - 2 * d)
' latitude arguments
S = F + (dL + 412 * SIN(2 * F) + 541 * SIN(LS)) / ARC
h = F - 2 * d
' latitude correction terms
N = -526 * SIN(h) + 44 * SIN(L + h) - 31 * SIN(h - L) - 23 * SIN(LS + h)
N = N + 11 * SIN(h - LS) - 25 * SIN(F - 2 * L) + 21 * SIN(F - L)
lmoon = p2 * fpart(L0 + dL / 1296000#) 'Lat in rads
bmoon = (18520# * SIN(S) + N) / ARC 'long in rads
' convert to equatorial coords using a fixed ecliptic
CB = COS(bmoon)
x = CB * COS(lmoon)
V = CB * SIN(lmoon)
W = SIN(bmoon)
y = COSEPS * V - SINEPS * W
Z = SINEPS * V + COSEPS * W
rho = SQR(1# - Z * Z)
dec = (360# / p2) * ATN(Z / rho)
ra = (48# / p2) * ATN(y / (x + rho))
IF ra < 0 THEN
ra = ra + 24#
END IF
END SUB

SUB quad (ym AS DOUBLE, y0 AS DOUBLE, yp AS DOUBLE, xe AS DOUBLE, ye AS DOUBLE, z1 AS DOUBLE, z2 AS DOUBLE, nz AS INTEGER)
' finds a parabola through three points and returns values of
' coordinates of extreme value (xe, ye) and zeros if any (z1, z2)
' assumes that the x values are -1, 0, +1
nz = 0
a = .5 * (ym + yp) - y0
b = .5 * (yp - ym)
c = y0
xe = -b / (2! * a) 'x coord of symmetry line
ye = (a * xe + b) * xe + c 'extreme value for y in interval
dis = b * b - 4! * a * c 'discriminant
IF dis > 0 THEN 'there are zeros
dx = .5 * SQR(dis) / ABS(a)
z1 = xe - dx
z2 = xe + dx
IF (ABS(z1) <= 1!) THEN nz = nz + 1 'This zero is in interval
IF (ABS(z2) <= 1!) THEN nz = nz + 1 'This zero is in interval
IF (z1 < -1!) THEN z1 = z2
END IF
END SUB

FUNCTION sinalt (iobj AS INTEGER, mjd0 AS DOUBLE, hour AS DOUBLE, glong AS DOUBLE, cphi AS DOUBLE, sphi AS DOUBLE)
' returns sine of the altitude of either the sun or the moon given the
' modified julian day number at midnight UT and the hour of the UT day,
' the longitude of the observer, and the sine and cosine of the latitude
' of the observer
ra = 0
dec = 0
instant = mjd0 + hour / 24#
t = (instant - 51544.5#) / 36525#
IF (iobj = 1) THEN
moon t, ra, dec
ELSE
sun t, ra, dec
END IF
tau = 15# * (lmst(instant, glong) - ra) 'hour angle of object
sinalt = sphi * sn(dec) + cphi * cn(dec) * cn(tau)
END FUNCTION

DEFSNG A-Z
FUNCTION sn# (x AS DOUBLE)
sn = SIN(x * .0174532925199433#)
END FUNCTION

DEFDBL A-Z
SUB sun (t AS DOUBLE, ra AS DOUBLE, dec AS DOUBLE)
' Returns RA and DEC of Sun to roughly 1 arcmin for few hundred
' years either side of J2000.0
p2 = 6.283185307#
COSEPS = .91748
SINEPS = .39778
m = p2 * fpart(.993133 + 99.997361# * t) 'Mean anomaly
dL = 6893# * SIN(m) + 72# * SIN(2 * m) 'Eq centre
L = p2 * fpart(.7859453# + m / p2 + (6191.2# * t + dL) / 1296000#)
' convert to RA and DEC - ecliptic latitude of Sun taken as zero
sl = SIN(L)
x = COS(L)
y = COSEPS * sl
Z = SINEPS * sl
rho = SQR(1# - Z * Z)
dec = (360# / p2) * ATN(Z / rho)
ra = (48# / p2) * ATN(y / (x + rho))
IF ra < 0 THEN ra = ra + 24
END SUB 

[DaveJJ]

@dodicat

I am p*$@ed off, ... you've looked at the code accepted all the $%!# in the statements and got it running. My hat is off!

The mjd thing mystified me because it looked like a global variable being used as a function label but it wasn't that obvious to me.

You've left the function declarations in place like function hm#(...) which i changed to function hm(........) as double. Does FB cope with this? and for instance, the .5 I have changed to 0.5 and stuff like .123456# I changed to 0.123456 to avoid problems with misinterpretation, was this a waste of time? If it was, FB is definitely a winner. I will try your suggestions.

many thanks

[DaveJJ]

@dodicat

I was struggling a bit last night so I left the code. This morning I tried to run the compiled code but got some errors referring to line numbers and function names which seemed strange.

Then I suddenly realised i was trying to run the QB text file and the compiler had created a file without the .BAS type appended, so all is ok. As an aside, I used fbc -lang qb ~/Freebasic/riset_mod.bas to compile. Is that ok or are there other suitable switches to use?

Thanks again

[dodicat]

#lang "qb"
#lang "fblite"
#lang "deprecated"

These can be written in your code, or used as a compiler switch.
E.G.
fbc -lang qb
as you have used.

Do you not use an IDE for FreeBasic?
E.G.
fbide.
Check out the help file for the three dialect options.
You also have options (among others) -gen gcc or -gen gas as a compiler switch,
The default is -gen gas which compiles via assembler, the -gen gcc compiles via gcc and the basic code is converted to C.
To use -gen gcc you should include the gcc executable download.

The long and short of it is, you have the options of coding in line number spaghetti code or traditional basic syntax on the one hand, or on the other hand, oop code optimized by Gcc, and a compiler slowly evolving towards C++ type power.

[DaveJJ]

Many thanks for that, I definitely need some practice.

The code had so many integer, double and single symbols strewn through it that I could not determine the type for new local variables and whether some had a global scope. Also, the code is difficult to decipher with no indentation at all. My modified code to turn it into something to work with FB does compile and run but does not return valid data. I guess it's easiest to keep it as QB but there are some very good additions to FB which I'd like to use.

I guess I'll raise a new thread about possibly mixing FB functions in with QB code.

Thanks again

[srvaldez]

you can find your program with indentation here http://www.stargazing.net/kepler/riset.bas

[DaveJJ]

A very useful version, thanks.

Whether it will help me is another question. It will be interesting to compare it with mine.