Go to Otherpower.com Home Page Go to Forcefield Shopping Cart Go to Wondermagnet.com Home Page
Front Page - [Homebrewed Electricity-- (wind) (solar) (hydro) (steam) (controls) (storage) (mechanical)] - Classifieds - Site News
Everything - Newbies - [Remote Living-- (housing) (heat) (light) (water)] - Rants & Opinion - Diaries - Our Products
software


By david williams, Section Rants & Opinion
Posted on Tue Jun 3rd, 2008 at 05:03:38 PM MST
proposal for a software collection

What would people think of having a section of the board devoted to relevant, public-domain (i.e. NOT copyright) computer software? I've seen lots of programs related to solar energy, wind-turbine design, Stirling engines, and other topics related to "other power". They are in various languages. Some use metric units, others not, and so on and so on. But they are all potentially useful. If we could make a collection that people could browse when they're looking for software, it might be widely appreciated.

Comments?

                      dow
software | 24 comments (24 topical, 0 editorial)

Re: software (3.00 / 0) (#1)
by Buzz Hacksaw (buzz_hacksaw@rogers.com) on Tue Jun 3rd, 2008 at 03:23:07 PM MST
(User Info)

I think it's a good idea. Go for it.
Buzz


Re: software (3.00 / 0) (#2)
by frackers (g8ecj at *nospam* gilks dot org) on Tue Jun 3rd, 2008 at 07:40:22 PM MST
(User Info) http://www.gilks.org

I wouldn't limit it to public domain - there are many license choices that allow an author to maintain copyright without restricting the free use of the software.

For example, all my software I release under the GPL Version 2. This would allow you to use it,  modify it if you wanted to (and distribute those changes if you wish) but if you put it into a proprietary product you have to either see me for a commercial license or publish your modifications so that the software remains 'free'. This allows my creative work to grow and be improved and allows me and others to benefit from the improvements others make.
Robin - Down Under (or are you Up Over)



Re: software (3.00 / 0) (#3)
by electrondady1 on Wed Jun 4th, 2008 at 10:08:58 AM MST
(User Info)

some were i saw a  a program that was designed to turn your old computor into am oscilloscope.
that would be slick.

[ Parent ]


Re: software (3.00 / 0) (#9)
by elvin1949 (elvin1949@yahoo.com) on Sat Jun 7th, 2008 at 12:05:00 AM MST
(User Info)

 I think that turns a small b/w tv into a scope.
Runs under win-xp.  I looked at it a few yr's ago.
 Not for me [ i don't use XP ]

 Google TV-oscilloscope   That should bring it up

Been a long time,hope this is right.

later
Elvin

[ Parent ]



Re: software (3.00 / 0) (#4)
by david williams on Thu Jun 5th, 2008 at 11:50:40 AM MST
(User Info)

Sounds good. I guess if an author uploads something on which he wants to retain copyright, he should say so clearly in the program or in text attached to it, and state what restrictions he wants to put on its distribution and use.

Any program that is uploaded by its author without such statements should be assumed to be public-domain, free for any kind of use.

Maybe, to be safe, we should allow uploads only by the author of the program. What do you think about that?

                            dow


[ Parent ]



Re: software (3.00 / 0) (#6)
by DamonHD (d@hd.org) on Thu Jun 5th, 2008 at 12:45:22 PM MST
(User Info) http://www.earth.org.uk/

The law in most places (eg all signatories to the Berne Convention) is quite clear that no announcement has to be made for copyright to remain intact, so assuming that something where you can't see the copyright notice is public-domain is very very dubious and is wishful thinking.

Just say No!

B^>

Rgds

Damon

[ Parent ]



Re: software (3.00 / 0) (#5)
by david williams on Thu Jun 5th, 2008 at 12:00:15 PM MST
(User Info)

Here's a program that should be useful in the context of solar energy. It does the astronomical and trigonometrical calculations that have to be done in a computer-controlled sun-tracker or heliostat. With the addition of code to allow the computer to control motors that turn the tracker or mirror, the program could become the software for such a machine.

I'll post two versions. The first is in QBasic but will also run under Borland Turbo BASIC and some other implementations. The second is in a very generic BASIC, with line numbers. It will run (and has been tested) under many varieties of BASIC. The QBasic version is the "real" program, which should be easier to read and understand than the other.

If anyone wants to translate it to another language and post the result, please go ahead!

This is public-domain, free for any kind of use, including commercial use.

Enjoy!

                               dow

------------------------------------------------

' SunAlign.BAS (Version for QBasic and similar dialects)

' Calculates position of sun in sky, as azimuth (compass bearing
' measured clockwise from True North) and angle of elevation, as
' seen from any place on earth, on any date and any time.
' Also calculates alignment of a heliostat mirror.

' David Williams
' P.O. Box 48512
' 3605 Lakeshore Blvd. West
' Toronto, Ontario. M8W 4Y6
' Canada

' Initially dated 2007 Jul 07
' This version 2008 Jan 13

' All angles in radians except in i/o routines DegIn and DegOut

DECLARE SUB C2P (X, Y, Z, AZ, EL)
DECLARE SUB P2C (AZ, EL, X, Y, Z)
DECLARE FUNCTION Ang (X, Y)
DECLARE SUB DegIn (P$, X)
DECLARE SUB DegOut (P$, X)

CONST PY = 3.1415926536# ' "PI" not assignable in some BASICs
CONST DR = 180 / PY ' degree / radian factor
W = 2 * PY / 365 ' earth's mean orbital angular speed in radians/day
WR = PY / 12' earth's speed of rotation relative to sun (radians/hour)
C = -23.45 / DR ' reverse angle of earth's axial tilt in radians
ST = SIN(C) ' sine of reverse tilt
CT = COS(C) ' cosine of reverse tilt
E2 = 2 * .0167 ' twice earth's orbital eccentricity
SN = 10 * W ' 10 days from December solstice to New Year (Jan 1)
SP = 12 * W ' 12 days from December solstice to perihelion

CLS

Menu:

PRINT "1. Calculate sun's position"
PRINT "2. Calculate mirror orientation"
PRINT "3. Calculate both"
PRINT "4. Quit program"
PRINT
PRINT "Which? (1 - 4)";

DO
  S% = VAL(INKEY$)
LOOP UNTIL S% >= 1 AND S% <= 4
PRINT S%

IF S% = 4 THEN END

' Note: For brevity, no error checks on user inputs

PRINT
PRINT "Use negative numbers for directions opposite to those shown."
PRINT
DegIn "Observer's latitude (degrees North)", LT
DegIn "Observer's longitude (degrees East)", LG
INPUT "Time Zone (+- hours from GMT/UT)"; TZN
INPUT "Time (HH,MM) (24-hr format)"; HR, MIN
INPUT "Date (M#,D#)"; Mth%, Day%

PRINT

CL = PY / 2 - LT ' co-latitude

D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365
' day of year (D = 0 on Jan 1)

A = W * D + SN ' orbit angle since solstice at mean speed
B = A + E2 * SIN(A - SP) ' angle with correction for eccentricity

C = (A - ATN(TAN(B) / CT)) / PY
SL = PY * (C - INT(C + .5))' solar longitude relative to mean position

C = ST * COS(B)
DC = ATN(C / SQR(1 - C * C)) ' solar declination (latitude)
' arcsine of C. ASN not directly available in QBasic

LD = (HR - TZN + MIN / 60) * WR + SL + LG ' longitude difference

CALL P2C(LD, DC, sX, sY, sZ) ' polar axis (perpend'r to azimuth plane)
CALL C2P(sY, sZ, sX, sAZ, sEL) ' horizontal axis
CALL P2C(sAZ - CL, sEL, sY, sZ, sX) ' rotate by co-latitude

IF sZ < 0 THEN
   BEEP
   PRINT "Sun Below Horizon"
   PRINT
   GOTO NewCalc
END IF

IF S% <> 2 THEN ' calculate and display sun's position

   CALL C2P(sX, sY, sZ, sAZ, sEL) ' vertical axis

   DegOut "Sun's azimuth: ", sAZ
   DegOut "Sun's elevation: ", sEL

   PRINT

END IF

IF S% > 1 THEN ' calculate and display mirror orientation

   PRINT "For target direction of light reflected from mirror:"
   DegIn "Azimuth of target direction (degrees)", tAZ
   DegIn "Elevation of target direction (degrees)", tEL

   PRINT

   CALL P2C(tAZ, tEL, tX, tY, tZ) ' target vector X,Y,Z
   CALL C2P(sX + tX, sY + tY, sZ + tZ, mAZ, mEL)
   ' angle bisection by vector addition

   PRINT "Mirror aim direction (perpendicular to surface):"
   DegOut "Azimuth: ", mAZ
   DegOut "Elevation: ", mEL

   PRINT

END IF

NewCalc:

PRINT
PRINT "New Calculation"
PRINT

GOTO Menu

FUNCTION Ang (X, Y)
 ' calculates angle from positive X axis to vector to (X,Y)

  SELECT CASE SGN(X)
    CASE 1: Ang = ATN(Y / X)
    CASE -1: Ang = ATN(Y / X) + PY
    CASE ELSE: Ang = SGN(Y) * PY / 2
  END SELECT

END FUNCTION

SUB C2P (X, Y, Z, AZ, EL)
  ' Cartesian to Polar. Convert from X,Y,Z to AZ,EL

  EL = Ang(SQR(X * X + Y * Y), Z)
  A = Ang(Y, X)
  IF A < PY THEN AZ = A + PY ELSE AZ = A - PY

END SUB

SUB DegIn (P$, X)
  ' Input angle in degrees and convert to radians

  PRINT P$;
  INPUT N
  X = N / DR

END SUB

SUB DegOut (P$, X)
 ' converts radians to degrees, rounds to nearest 0.1, and prints

  S$ = LTRIM$(STR$(INT(10 * ABS(X * DR) + .5)))
  IF S$ = "3600" THEN S$ = "0"
  IF LEN(S$) = 1 THEN S$ = "0" + S$
  IF X < 0 THEN IF VAL(S$) THEN S$ = "-" + S$
  PRINT P$; LEFT$(S$, LEN(S$) - 1); "."; RIGHT$(S$, 1); " degrees"

END SUB

SUB P2C (AZ, EL, X, Y, Z)
  ' Polar to Cartesian. Convert from AZ,EL to X,Y,Z

  Z = SIN(EL)
  C = -COS(EL)
  X = C * SIN(AZ)
  Y = C * COS(AZ)

END SUB

--------------------------------------------------------------

--------------------------------------------------------------

  1. REM SunAlign.BAS (Generic BASIC version)
  2. REM Calculates position of sun in sky, as azimuth (compass bearing
  3. REM measured clockwise from True North) and angle of elevation, as
  4. REM seen from any place on earth, on any date and any time.
  5. REM Also calculates alignment of a heliostat mirror.
  6. REM David Williams
  7. REM P.O. Box 48512
  8. REM 3605 Lakeshore Blvd. West
  9. REM Toronto, Ontario. M8W 4Y6
  10. REM Canada
  11. REM Original date 2007 Jul 07. This version 2007 Oct 07
  12. REM Note: For brevity, no error checks on user inputs
  13. CLS
  14. PRINT "Use negative numbers for opposite directions."
  15. INPUT "Observer's latitude (degrees North)"; LT
  16. INPUT "Observer's longitude (degrees East)"; LG
  17. INPUT "Date (M#,D#)"; Mth, Day
  18. INPUT "Time (HH,MM) (24-hr format)"; HR, MIN
  19. INPUT "Time Zone (+- hours from GMT/UT)"; TZN
  20. PY = 4 * ATN(1): REM "PI" not assignable in some BASICs
  21. DR = 180 / PY: REM degree/radian factor
  22. W = 2 * PY / 365: REM earth's mean orbital speed in radians/day
  23. C = -23.45 / DR: REM reverse angle of axial tilt in radians
  24. ST = SIN(C): REM sine of reverse tilt
  25. CT = COS(C): REM cosine of reverse tilt
  26. E2 = 2 * .0167: REM twice earth's orbital eccentricity
  27. SP = 12 * W: REM 12 days from December solstice to perihelion
  28. D = INT(30.6 * ((Mth + 9) MOD 12) + 58.5 + Day) MOD 365
  29. A = W * (D + 10): REM Solstice 10 days before Jan 1
  30. B = A + E2 * SIN(A - SP)
  31. C = (A - ATN(TAN(B) / CT)) / PY
  32. ET = 720 * (C - INT(C + .5)): REM equation of time
  33. REM in 720 minutes, earth rotates PI radians relative to sun
  34. C = ST * COS(B)
  35. EL = ATN(C / SQR(1 - C * C)) * DR: REM solar declination
  36. AZ = 15 * (HR - TZN) + (MIN + ET) / 4 + LG: REM longitude diff
  37. GOSUB 800
  38. R = SQR(Y * Y + Z * Z)
  39. AX = Y: AY = Z: GOSUB 710
  40. A = AA + (90 - LT) / DR
  41. Y = R * COS(A)
  42. Z = R * SIN(A)
  43. GOSUB 740
  44. PRINT : REM AZ & EL are now sun's azimuth & elevation in degrees
  45. IF EL < 0 THEN PRINT "Sun Below Horizon": END
  46. R = AZ: GOSUB 870: PRINT "Sun's azimuth: "; R; " degrees"
  47. R = EL: GOSUB 870: PRINT "Sun's elevation: "; R; " degrees"
  48. PRINT
  49. INPUT "Calculate heliostat mirror alignment (y/n)"; K$
  50. IF K$ = "N" OR K$ = "n" THEN END
  51. SX = X: SY = Y: SZ = Z
  52. PRINT
  53. INPUT "Azimuth of target direction (degrees)"; AZ
  54. INPUT "Elevation of target direction (degrees)"; EL
  55. GOSUB 800
  56. X = X + SX: Y = Y + SY: Z = Z + SZ: GOSUB 740
  57. PRINT : REM AZ & EL are now aim azimuth & elevation in degrees
  58. PRINT "Mirror aim direction (perpendicular to surface):"
  59. R = AZ: GOSUB 870: PRINT "Azimuth: "; R; " degrees"
  60. R = EL: GOSUB 870: PRINT "Elevation: "; R; " degrees"
  61. END
  62. IF AX = 0 THEN AA = SGN(AY) * PY / 2: RETURN
  63. AA = ATN(AY / AX): IF AX < 0 THEN AA = AA + PY
  64. RETURN
  65. AX = SQR(X * X + Y * Y): AY = Z: GOSUB 710
  66. EL = AA * DR
  67. AX = Y: AY = X: GOSUB 710
  68. AZ = AA * DR
  69. IF AZ < 180 THEN AZ = AZ + 180 ELSE AZ = AZ - 180
  70. RETURN
  71. E = EL / DR
  72. A = AZ / DR
  73. Z = SIN(E)
  74. C = 0 - COS(E): REM Won't work without "0" in Liberty Basic
  75. X = C * SIN(A)
  76. Y = C * COS(A)
  77. RETURN
  78. R = INT(10 * R + .5): IF R = 3600 THEN R = 0
  79. R = R / 10
  80. RETURN
--------------------------------------------------



Re: software (3.00 / 0) (#7)
by TomW on Thu Jun 5th, 2008 at 01:29:36 PM MST
(User Info)

Yeah, anyone who can translate this to a shell or perl script or other unix tool[s] that would allow me to locate the solar position based on my day/time and timezone would be cool beans. I can track horizon to horizon [treetops really] with a resolution of 515 steps. Roughly 3 degree steps.

Tom

PS

Speaking of solar, I just read and measured 56 inches multiple times then promptly laid it out as 59 while upgrading my trackers rack. Bolts are at least easy to change.

"Education consists mainly of what we have unlearned."--Mark Twain
[ Parent ]



Re: software (3.00 / 0) (#10)
by david williams on Sat Jun 7th, 2008 at 11:08:49 AM MST
(User Info)

If you, or anyone else, want to translate my Sunalign program into a perl script or anything else, I'd be happy to help however I can. I do undertand the calculation. I could do it with a pocket calculator if I had to. I coded it into BASIC because that's the only programming language I'm any good at. If anyone wants to re-code it into any other language (providing, of course, that it has provision for floating-point math and trig functions), and doesn't understand it from the BASIC, just contact me and I'll try to explain how it works. Contact me here or, probably better, send me an e-mail:

david.williams@bayman.org

Later...
                              dow


[ Parent ]



Re: software (3.00 / 0) (#11)
by DamonHD (d@hd.org) on Sat Jun 7th, 2008 at 11:20:33 AM MST
(User Info) http://www.earth.org.uk/

It's potentially the sort of thing that might be coded up as JavaScript on pages hosted here at FL/OP...

Easy free service which costs FL nothing (except in TomW et al handling newbie impotent rage when it doesn't give the desired answer...)

Rgds

Damon

[ Parent ]



Re: software (3.00 / 0) (#12)
by david williams on Tue Jun 10th, 2008 at 12:16:51 PM MST
(User Info)

Translating Sunalign to JavaScript sounds like a good idea. Making it so it doesn't work doesn't! There are already other programs out there that are claimed to do the same thing, but they aren't always reliable. Some, for example, don't work correctly for locations in the southern hemisphere. I am very confident that Sunalign always works correctly, with accuracy that is more than sufficient for solar-energy purposes. Let's keep it that way. Please don't post any version that doesn't work reliably.

                                dow


[ Parent ]



Re: software (3.00 / 0) (#13)
by DamonHD (d@hd.org) on Tue Jun 10th, 2008 at 02:28:12 PM MST
(User Info) http://www.earth.org.uk/

Your point is well made, though I was referring in part to the kind of rage incurred when the program gives the mathematically correct answer, but not the answer hoped for by the innumerate and physics-detesting user who wanted to stick it to the (power) man...

Rgds

Damon

[ Parent ]



Re: software (3.00 / 0) (#14)
by david williams on Wed Jun 11th, 2008 at 02:56:11 PM MST
(User Info)

In the middle of a Canadian winter, I'd love the sun to be high in the sky. I could easily write a program to tell me that this is so, but I don't think it would help much!

                            dow


[ Parent ]



Re: software (3.00 / 0) (#15)
by david williams on Fri Aug 15th, 2008 at 02:24:55 PM MST
(User Info)

I know of someone who has translated my BASIC Sunalign program into Perl. He says it seems to work fine, but he wants to give it a very thorough testing and polishing before releasing it into the big wide world.

When he feels it's ready, I'll ask him to post it here.


[ Parent ]



Re: software (3.00 / 0) (#16)
by bills on Wed Aug 20th, 2008 at 07:20:01 AM MST
(User Info)

Tom,

I have ported David Williams' SunAlign to a Perl program although I haven't completed the comprehensive testing on it yet. If you would like a copy drop me an email and I will send it on. Maybe you could try it and let me know what you think? I use gmail. [First initial][middle initial][last name]

William K Stewart (Bill)

[ Parent ]



Re: software (3.00 / 0) (#8)
by david williams on Fri Jun 6th, 2008 at 01:44:28 PM MST
(User Info)

Here's a solar "metaprogram". It's really a program to aid in writing programs for solar-energy applications, and also for other solar stuff such as sundials. It contains a subroutine, a QBasic FUNCTION, that calculates the sun's declination (latitude) and the Equation of Time on any date. The Equation of Time is the difference between the time shown on a sundial and that shown by a clock. It is essentially the sun's longitude, relative to its mean position. One degree of longitude equals four minutes in the Equation of Time.

As I said, it's in QBasic. It contains some graphics, so it can't easily be completely translated to another dialect, but the FUNCTION coding should be easy to translate.

It's public-domain. Free for any kind of use, including commercial use.

I hope someone finds it useful.

                             dow

---------------------------------------------------------

' ETIMSDEC.BAS
' David Williams, 2003

' P.O. Box 48512
' 3605 Lakeshore Blvd. West
' Toronto, Ontario. M8W 4Y6
' Canada

' This version dated 2007 Jun 13

DECLARE FUNCTION ET.Dec (D, F%)
DECLARE FUNCTION ROff$ (X)
DECLARE SUB Waitkey (X%)

CL0 = 7: CL1 = 12: CL2 = 15   ' colours used

SCREEN 12

COLOR CL2
CLS

PRINT "ETIMSDEC.BAS"
PRINT
PRINT "Shows Equation of Time and Solar Declination calculations"
PRINT "(performed in function ET.Dec) and compares their results"
PRINT "graphically with published values, showing close agreement."
PRINT
PRINT "These calculations can be used in the programs of computer-"
PRINT "controlled solar-energy equipment, such as sun trackers and"
PRINT "heliostats."
PRINT
PRINT "Equation of Time is difference between Solar Time and Mean"
PRINT "Time. Sundials show solar time. Clocks show mean time. The"
PRINT "Equation of Time is related to solar longitude, relative to"
PRINT "its mean value (the value it would have if the sun's apparent"
PRINT "motion were at uniform speed). One degree westward of solar"
PRINT "longitude corresponds to four minutes in positive direction in"
PRINT "the Equation of Time."
PRINT
PRINT "Solar Declination is latitude of sun in celestial coordinates."
PRINT

' initialize array used in graphing routine
DIM ML(1 TO 13)
FOR X = 1 TO 13
  READ ML(X)
NEXT
DATA 33,29,33,31,33,31,33,33,31,33,31,33,0
' Month lengths adjusted for numbers of pixels on screen

DO

  COLOR CL0

  PRINT "1. Find Equation of Time and Solar Declination on given date"
  PRINT "2. Draw Equation of Time Graph"
  PRINT "3. Draw Solar Declination Graph"
  PRINT "4. Draw Solar Analemma"
  PRINT "5. Quit program"
  PRINT
  PRINT "Which (1 - 5)? ";
  WHILE INKEY$ <> "": WEND
  DO
    K$ = INKEY$
  LOOP UNTIL K$ >= "1" AND K$ <= "5"
  PRINT K$

  SELECT CASE VAL(K$)
    CASE 1: GOSUB Calc
    CASE 2, 3: GOSUB Graph
    CASE 4: GOSUB Analemma
    CASE 5: EXIT DO
  END SELECT

LOOP

GOTO WayOut

Calc:

PRINT
COLOR CL1

DO
  INPUT "Date (M#,D#)"; Mth%, Day%
LOOP UNTIL Mth% > 0 AND Mth% < 13 AND Day% > 0 AND Day% < 32

D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365
' day of 365-day year. Jan 1 = 0. Dec 31 = 364

PRINT
PRINT "Day number"; D + 1; "of year" ' conventional range, 1 - 365

ET = ET.Dec(D, 1) ' Equation of Time
DC = ET.Dec(D, 0) ' Declination

PRINT
PRINT "Equation of Time: "; ROff$(ET); " minutes"
PRINT "Solar Declination: "; ROff$(DC); " degrees"
PRINT

RETURN ' to menu

Graph:

F% = (K$ = "2") ' true if Equation of Time selected

' Set up graph
CLS

IF F% THEN
  PRINT TAB(30); "EQUATION OF TIME"
  LOCATE 3, 15
  PRINT "Graph shows difference in minutes between clock"
  PRINT TAB(15); "and sundial time. Positive difference means"
  PRINT TAB(15); "sundial is ahead of clock, and vice versa."
  LOCATE 24, 15
  COLOR CL1
  PRINT "Graph shows results of calculation."
  COLOR CL2
  PRINT TAB(15); "Crosses show published values of Equation of Time."
  RESTORE ETData
ELSE
  PRINT TAB(25); "Declination of Sun, in Degrees"
  LOCATE 3, 10
  COLOR CL1
  PRINT "Graph shows calculated function. ";
  COLOR CL2
  PRINT "Crosses show published values."
  RESTORE DeclnData
END IF

READ ML%, VT%, VB%, LL%, UL%

COLOR CL0
LOCATE 16, 67
PRINT "-="

FOR T = LL% TO UL% STEP 5
  LINE (137, 247 - 6.4 * T)-(530, 247 - 6.4 * T)
  LOCATE 16 - T / 2.5, 14
  PRINT RIGHT$(" " + STR$(T), 3);
  IF T = 0 THEN PRINT " ="
NEXT
LOCATE ML%, 20
PRINT "Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec"

' Do graphing
T = 0
D = 0
MC = 0
FOR M = 0 TO 384
  M1 = 147 + M
  IF M = T THEN
    ' Draw vertical line
    COLOR CL0
    LINE (M1, VT%)-(M1, VB%)
    MC = MC + 1
    T = T + ML(MC)
    MS = 0
  ELSE
    ' Plot point(s)
    COLOR CL1
    YY = 247 - 6.4 * ET.Dec(D, F%) ' function used
    IF POINT(M1, YY) <> CL2 THEN PSET (M1, YY) ' don't overwrite cross
    IF MS = 0 OR MS = 10 OR MS = 20 THEN
      ' Plot cross representing published value
      COLOR CL2
      READ E ' published value
      E = FIX(E) + (E - FIX(E)) * 5 / 3 ' see note after E of T data
      YY = 247 - 6.4 * E
      LINE (M1 - 2, YY)-(M1 + 2, YY) ' lines make
      LINE (M1, YY - 2)-(M1, YY + 2) ' cross
    END IF
    MS = MS + 1
    D = D + 1
  END IF
  IF ML(MC) = 33 AND M = T - 18 THEN M = M + 1
NEXT

Waitkey 25

RETURN ' to menu

ETData:

DATA 7, 97, 343, -15, 20: ' plotting parameters for E. of T. graph

DATA -3.12,-7.38,-11.08,-13.33,-14.19,-13.49,-12.34,-10.18,-7.28
DATA -4.08,-1.16,1,2.51,3.4,3.34,2.25,.39,-1.28
DATA -3.33,-5.16,-6.15,-6.16,-5.14,-3.16,-.12,3.08,6.4
DATA 10.05,13.02,15.12,16.2,16,14.16,11.11,7.02,2.13
' Above data show published values of Equation of Time on
' 1st, 11th and 21st days of months. Decimal point separates
' minutes, to left, from seconds, to right, e.g. 3.16 means
' 3 minutes and 16 seconds. Seconds part is multiplied by 5/3
' to get fractional minutes after READ instruction.

DeclnData:

DATA 5, 65, 407, -25, 25: ' plotting parameters for Decl'n graph

DATA -23.04,-21.56,-20.05,-17.2,-14.18,-10.52,-7.49,-3.57,0
DATA 4.18,8.07,11.39,14.54,17.43,20.04,21.58,23.02,23.26
DATA 23.09,22.11,20.36,18.1,15.27,12.19,8.3,4.48,.57
DATA -2.57,-6.48,-10.29,-14.14,-17.15,-19.47,-21.43,-22.57,-23.26
' Above data show published values (in degrees and minutes) of solar
' declination on 1st, 11th and 21st days of month. See note after
' Equation of Time data.

' Published values are from the book:
' Sundials, Their Theory and Construction
' by: Albert E. Waugh
' Dover Publications, New York, 1973
' ISBN 0-486-22947-5

Analemma:

CLS

RESTORE Anadata

PRINT "SOLAR ANALEMMA"
PRINT
PRINT
PRINT "Shows apparent path of sun,"
PRINT "through the year, relative"
PRINT "to its mean position"
PRINT "(marked by cross)."

LOCATE 3, 60
PRINT "N"
LOCATE 7, 60
PRINT "S"
LOCATE 5, 56
PRINT "E       W"

LINE (475, 52)-(475, 90)  ' compass lines
LINE (456, 71)-(494, 71)

LOCATE 9, 1
PRINT TAB(53); "Looking at sky,"
PRINT TAB(53); "if North is upward,"
PRINT TAB(53); "East is to left."

LOCATE 24, 1
PRINT "For dimensions, see graphs"
PRINT "of Solar Declination and"
PRINT "Equation of Time."
PRINT "+4 minutes in Equation of Time"
PRINT "equal 1 degree Westward."

LOCATE 15, 1
PRINT TAB(50); "Celestial"
PRINT TAB(50); "Equator"

X0 = 320  ' origin coordinates
Y0 = 240

LINE (X0, Y0 + 5)-(X0, Y0 - 5) ' cross at (0,0)
LINE (X0 + 5, Y0)-(X0 - 5, Y0)

FOR S = -1 TO 1 STEP 2
  LINE (X0, Y0 + S * 210)-(X0, Y0 + S * 225) ' lines marking N-S
  LINE (X0 + S * 30, Y0)-(X0 + S * 55, Y0)   ' and E-W
NEXT

LOCATE 10, 1
PRINT "Month colours:";

M = 0 ' month number
DM = 0 ' day of month
LM = 0 ' length of month

' draw analemma
FOR DY = 0 TO 364 ' day of year
  IF DM = LM THEN ' start new month
    M = M + 1
    C = M XOR 14 * (M AND 1) ' twiddle bits to get good contrasts
    IF C = 8 THEN C = 14  ' yellow better than dull grey
    COLOR C
    READ LM, MN$
    PRINT TAB(17); MN$ ' month name in correct colour
    DM = 0 ' reset day of month
  END IF
  DM = DM + 1
  FOR D = DY TO DY + .9 STEP .25 ' 4 points per day to avoid gaps
    E = ET.Dec(D, 1) ' equation of time
    L = ET.Dec(D, 0) ' solar declination (latitude)
    PSET (X0 + 2.5 * E, Y0 - 10 * L) ' 4:1 ratio.  4 mins = 1 degree
  NEXT
NEXT

Waitkey 48

RETURN ' to menu

Anadata:

DATA 31,Jan,28,Feb,31,Mar,30,Apr,31,May,30,Jun
DATA 31,Jul,31,Aug,30,Sep,31,Oct,30,Nov,31,Dec
' month lengths and names

WayOut:

COLOR CL2
CLS
PRINT "These calculations of Equation of Time and Solar Declination"
PRINT "are simplified and approximate. However, they are quite good."
PRINT "All differences between calculated and published values of"
PRINT "the solar declination are small compared with the angular"
PRINT "size of the sun in the sky. Similarly, differences between"
PRINT "calculated and published values of the equation of time are"
PRINT "small compared with the time the sun takes to traverse its"
PRINT "own diameter as it moves across the sky. The non-zero size"
PRINT "of the sun, rather than any inaccuracies of calculation,"
PRINT "is the limiting factor on how accurately the sun can be"
PRINT "tracked in solar energy applications, using this routine."
PRINT
PRINT "Leap years cause minor fluctuations in the Equation of Time"
PRINT "and Solar Declination on given dates of the year. However,"
PRINT "the errors that are introduced by ignoring these fluctuations"
PRINT "are always smaller than the apparent radius of the sun. This"
PRINT "program therefore ignores leap years. February 29, when it"
PRINT "occurs, is considered to be the same day as March 1."
PRINT
PRINT "The Equation of Time is treated here as the correction to be"
PRINT "subtracted from a sundial reading to get local mean (clock)"
PRINT "time. It is therefore positive when sundial is ahead of clock."
PRINT "This is the usual sign convention, but the opposite usage is"
PRINT "sometimes found. Take care when comparing values of the"
PRINT "Equation of Time from different sources."

Waitkey 1

PRINT "Note that atmospheric refraction of light affects the apparent"
PRINT "position of the sun in the sky when it is close to the"
PRINT "horizon. Sunlight is then too weak to be used for most"
PRINT "solar-energy applications, so refraction is not usually"
PRINT "included in calculations related to solar energy. This program"
PRINT "does not take account of refraction."
PRINT
PRINT "Comments in the function ET.Dec show how quantities in it are"
PRINT "related to astronomical parameters. Over long periods of time"
PRINT "(many decades or more), some of these parameters will change"
PRINT "significantly. The function should be modified accordingly if"
PRINT "it is to be used over very long periods of time."

Waitkey 1

END

FUNCTION ET.Dec (D, F%) STATIC
  ' Calculates equation of time, in minutes, or solar declination,
  ' in degrees, on day number D of year. (D = 0 on January 1.)
  ' F% selects function: True (non-zero) for Equation of Time,
  ' False (zero) for Declination.
  ' STATIC means variables are preserved between calls of function

IF PI = 0 THEN ' first call, initialize constants

 PI = 4 * ATN(1)
 W = 2 * PI / 365 ' earth's mean orbital angular speed in radians/day
 DR = 180 / PI ' degree/radian factor
 C = -23.45 / DR ' reverse angle of earth's axial tilt in radians
 ST = SIN(C) ' sine of reverse tilt
 CT = COS(C) ' cosine of reverse tilt
 E2 = 2 * .0167 ' twice earth's orbital eccentricity
 SP = 12 * W ' 12 days from December solstice to perihelion
 D1 = -1 ' holds last D. Saves time if D repeated for both functions

END IF

IF D <> D1 THEN ' new value of D
  A = W * (D + 10) ' Solstice 10 days before Jan 1
  B = A + E2 * SIN(A - SP)
  D1 = D
END IF

IF F% THEN  ' equation of time calculation
  C = (A - ATN(TAN(B) / CT)) / PI
  ET.Dec = 720 * (C - INT(C + .5))
  ' in 720 minutes, earth rotates PI radians relative to sun

ELSE ' declination calculation
  C = ST * COS(B)
  ET.Dec = ATN(C / SQR(1 - C * C)) * DR
  ' arcsine of C in degrees. ASN not directly available in QBasic

END IF

END FUNCTION

FUNCTION ROff$ (X)
 ' neatly rounds number to one place of decimals

S$ = LTRIM$(STR$(INT(10 * ABS(X) + .5)))
IF LEN(S$) = 1 THEN S$ = "0" + S$
IF VAL(S$) THEN S$ = MID$("-x+", SGN(X) + 2, 1) + S$
ROff$ = LEFT$(S$, LEN(S$) - 1) + "." + RIGHT$(S$, 1)

END FUNCTION

SUB Waitkey (X%)
 ' prints prompt at position X% on bottom row, then waits

COLOR 15
LOCATE 30, X%
PRINT "* Press any key to continue *";
WHILE INKEY$ <> "": WEND
WHILE INKEY$ = "": WEND
CLS

END SUB

---------------------------------------------------



Re: software (3.00 / 0) (#18)
by david williams on Tue Aug 26th, 2008 at 01:31:21 PM MST
(User Info)

Thanks Bill. It looks good. Strange, to my eyes, but good.


Re: software (3.00 / 0) (#19)
by bills on Wed Sep 3rd, 2008 at 07:20:50 AM MST
(User Info)

Thanks David, It looks like a bunch of the symbols have disappeared from the code in the prints there should be backslash n for a newline when I posted it they seem to have been stripped. I'll have to look into this.
I will try to re-post a corrected copy.

Bill

[ Parent ]



Re: software (3.00 / 0) (#20)
by bills on Wed Sep 3rd, 2008 at 07:26:56 AM MST
(User Info)

#!/usr/bin/perl

use strict;
use warnings;

use Math::Trig;

#
# SunAlign.pl (Version for Perl)
#
# Calculates position of sun in sky, as azimuth (compass bearing
# measured clockwise from True North) and angle of elevation, as
# seen from any place on earth, on any date and any time.
# Also calculates alignment of a heliostat mirror.
#
# David Williams
# P.O. Box 48512
# 3605 Lakeshore Blvd. West
# Toronto, Ontario. M8W 4Y6
# Canada
#
# Initially dated 2007 Jul 07 (Basic version)
# This version 2008 August 08
#
# Perl version ported by William K Stewart, 2008 Aug 15
# gmail username = lc([first initial][middle initial][lastname])
#
# All angles in radians except in i/o routines DegIn and DegOut
#

my $PY = 3.1415926536;           # "PI" not assignable in some BASICs
my $DR = 180/$PY;                # degree / radian factor
my $W  = 2 * $PY/365;            # earth's mean orbital angular speed in radians/day
my $WR = $PY/12;                 # earth's speed of rotation relative to sun (radians/hour)
my $C  = -23.45/$DR;             # reverse angle of earth's axial tilt in radians
my $ST = sin($C);                # sine of reverse tilt
my $CT = cos($C);                # cosine of reverse tilt
my $E2 = 2 * 0.0167;             # twice earth's orbital eccentricity
my $SN = 10 * $W;                # 10 days from December solstice to New Year (Jan 1)
my $SP = 12 * $W;                # 12 days from December solstice to perihelion

my($LT, $LG, $tAZ, $tEL) = ();

while(1) {
   my $S;
    do{
        print "1. Calculate sun's position\n";
        print "2. Calculate mirror orientation\n";
        print "3. Calculate both\n";
        print "4. Quit program\n\n";
        print "Which? (1 - 4): ";
        $S = <STDIN>;
        chomp $S;
    } while(!($S >= 1 && $S <= 4));
    print "\nYou selected: $S\n";

    last if $S == 4;

    print "\nUse negative numbers for directions opposite to those shown.\n\n";

    DegIn ("Observer's latitude (degrees North)", \$LT);
    DegIn ("Observer's longitude (degrees East)", \$LG);

    my $TZN;
    do {
        print "Time Zone (+- hours from GMT/UT): ";
        $TZN = <STDIN>;
        chomp $TZN;
    } while($TZN !~ /^-?\d+$/);

    my $time;
    do {
        print "Time (HH,MM) (24-hr format): ";
        $time = <STDIN>;
        chomp $time;
    } while($time !~ /^\d+,\d+$/);
    my($HR, $MIN) = split(",", $time);

    my $date;
    do {
        print "Date (M#,D#): ";
        $date = <STDIN>;
        chomp $date;
    } while($date !~ /^\d+,\d+$/);
    my($Mth, $Day) = split(",", $date);

    print "\n";
 # co-latitude
    my $CL = $PY/2 - $LT;

    #day of year (D = 0 on Jan 1)
    my $D = (30.6 * (($Mth + 9) % 12) + 58.5 + $Day) % 365;

    # orbit angle since solstice at mean speed
    my $A = $W * $D + $SN;

    # angle with correction for eccentricity
    my $B = $A + $E2 * sin($A - $SP);

    $C = ($A - atan2(tan($B),$CT))/$PY;

    my $tmp = $C + .5;
    my $tmp1 = $C - 1;

    # solar longitude relative to mean position
    #my $SL = $PY * ($C - INT($C + .5));
    my $SL = $PY * ($C - int($C + .5));

    $C = $ST * cos($B);

    # solar declination (latitude)
    my $DC = atan2($C, sqrt(1 - $C * $C));

    # longitude difference
    my $LD = ($HR - $TZN + $MIN / 60) * $WR + $SL + $LG;

    # polar axis (perpend'r to azimuth plane)
    my($sX, $sY, $sZ);
    P2C($LD, $DC, \$sX, \$sY, \$sZ);

    # horizontal axis
    my($sAZ, $sEL);
    C2P($sY, $sZ, $sX, \$sAZ, \$sEL);

    # rotate by co-latitude
    P2C($sAZ - $CL, $sEL, \$sY, \$sZ, \$sX);

    if($sZ < 0) {
       print "\n\nSun Below Horizon\n\n";
       print "New Calculation\n";
       next;
    }
 if($S != 2) {
        # calculate and display sun's position
        # vertical axis
        C2P($sX, $sY, $sZ, \$sAZ, \$sEL);

        DegOut("Sun's azimuth: ", \$sAZ);
        DegOut("Sun's elevation: ", \$sEL);
        print "\n";
    }

    if($S > 1) {
        # calculate and display mirror orientation

        print "For target direction of light reflected from mirror:\n";
        DegIn ("Azimuth of target direction (degrees)", \$tAZ);
        DegIn ("Elevation of target direction (degrees)", \$tEL);

        print "\n";

        my($tX, $tY, $tZ);
        # target vector X,Y,Z
        P2C($tAZ, $tEL, \$tX, \$tY, \$tZ);

        my($mAZ, $mEL);
        # angle bisection by vector addition
        C2P($sX + $tX, $sY + $tY, $sZ + $tZ, \$mAZ, \$mEL);

        print "Mirror aim direction (perpendicular to surface):\n";
        DegOut("Azimuth: ", \$mAZ);
        DegOut("Elevation: ", \$mEL);
        print "\n";
    }
}

sub Ang {
    # calculates angle from positive X axis to vector to (X,Y)
    my($X, $Y) = @_;

    my $angle;
    my $val = sgn($X);

    if($val == 1) {
        $angle = atan($Y/$X);
    }
    elsif($val == -1) {
        $angle = atan($Y/$X) + $PY;
    }
    else {
        $angle = sgn($Y) * $PY/2;
    }

    return $angle;
}

sub C2P {
    # Cartesian to Polar. Convert from X,Y,Z to AZ,EL
    my($X, $Y, $Z, $AZ, $EL) = @_;

    $$EL = Ang(sqrt($X * $X + $Y * $Y), $Z);

    my $A = Ang($Y, $X);

    if($A < $PY) {
        $$AZ = $A + $PY;
    }
    else {
        $$AZ = $A - $PY;
    }
    return;
}

sub DegIn {
    # Input angle in degrees and convert to radians
    my($P, $X) = @_;

    my $N;
    do {
        print $P,": ";
        $N = <STDIN>;
        chomp $N;
    } while($N !~ /^-?\d+$/);

    $$X = $N/$DR;
    return;
}

sub DegOut {
    # converts radians to degrees, rounds to nearest 0.1, and prints
    my($P, $X) = @_;

    #my $S = LTRIM(STR(INT(10 * abs($$X * $DR) + .5)));
    my $S = 10 * abs($$X * $DR);
    $S = round($S);

    $S = 0 if($S eq "3600");
   $S = "0" . $S if(length($S) == 1);

    if($$X < 0) {
        if($S) {
            $S = "-" . $S;
        }
    }
    $S =~ s/(.+)(.)$/$1\.$2/;
    print "$P $S degrees\n";
    return;
}

sub P2C {
    # Polar to Cartesian. Convert from AZ,EL to X,Y,Z
    my($AZ, $EL, $X, $Y, $Z) = @_;

    $$Z = sin($EL);
    $$X = -cos($EL) * sin($AZ);
    $$Y = -cos($EL) * cos($AZ);
    return;
}

sub sgn {
    my($x) = @_;

    if($x < 0) {
        return -1
    }
    elsif($x == 0) {
        return 0;
    }
    else {
        return 1;
    }
}

sub round {
    my($number) = @_;
    return int($number + .5 * ($number <=> 0));
}

[ Parent ]



Re: software (3.00 / 0) (#21)
by david williams on Wed Sep 3rd, 2008 at 11:36:50 AM MST
(User Info)

Yes. It can be tricky to paste program code into a text-oriented system like this one. The system tries to make it look pretty, which may spoil the functioning of the program. Also, things may be wrongly interpreted as commands. For example, in SunAlign I have a variable called TZN (time zone). Some systems take TZN to be a command to switch to italic font, so weird things happen when the program is listed.

If anyone wants a "fresh" copy of the BASIC code for SunAlign, I'll e-mail it to them. Just write to me at:

david.williams@bayman.org

I'm sure you would do the same.

David


[ Parent ]



Re: software (3.00 / 0) (#22)
by TomW on Wed Sep 3rd, 2008 at 12:55:46 PM MST
(User Info)

david;

I got the last version posted by bills to run on the MacBook I couldn't get the first perl script post to work anywhere but didn't try to debug it far.

Thanks all around.

Tom

"Education consists mainly of what we have unlearned."--Mark Twain
[ Parent ]



Re: software (3.00 / 0) (#24)
by bills on Fri Sep 26th, 2008 at 12:40:41 PM MST
(User Info)

Tom, any others, the first version I posted, I posted incorrectly and the backslashes and some other symbols got filtered out. Please only use the second Perl version posted.

Best,
Bill

If the admin wanted to delete the first Perl version I posted it would probably eliminate future confusion.

[ Parent ]



Re: software (3.00 / 0) (#25)
by TomW on Fri Sep 26th, 2008 at 01:54:50 PM MST
(User Info)

bills;

Thanks for posting the working version. I removed the mangled version per your request.

Tom

"Education consists mainly of what we have unlearned."--Mark Twain
[ Parent ]



Re: software (3.00 / 0) (#23)
by david williams on Mon Sep 8th, 2008 at 11:54:59 AM MST
(User Info)

I took a hard look at the BASIC versions of SunAlign above. The QBasic version has weird italic in it, but that may not stop it from working. The generic-BASIC version has somehow acquired a period after every line number. I don't know why the system decided to do that. Unless the periods are removed, the code probably won't work. If you'd like me to send you a copy without any corruptions, just send me an e-mail.

dow

david.williams@bayman.org


[ Parent ]



software | 24 comments (24 topical, 0 editorial)
Display: Sort:
Menu
· create account
· How to use the board
· FAQs
· search the board
· Google search the board
· Old Otherpower Board

Login
Make a new account
Username:
Password:

Total Views
  113 Scoop users have viewed this posting.

Related Links
· Also by david williams

Powered by Scoop
You must be a registered user to post here. It's easy and free, and the link is on the upper right side of your page.
All trademarks and copyrights on this page are owned by their respective companies. Postings are owned by the poster, but may be deleted or moved at the ADMIN's sole discretion. The Rest © 2003 Forcefield.
You can Email the board ADMIN here. PLEASE include the username you signed up with!