'*******************************************************************
'*   Hor_analem #3, a DeltaCad macro for producing a horizontal    *
'*         analemmatic sundial with Time Zone Correction           * 
'* Copyright 2000 Fer de Vries & The North American Sundial Society*
'*                     http://sundials.org                         *
'*  This program may be circulated and modified as long as this    *
'*                      header remains intact.                     *
'*******************************************************************


'Formulae
'Choose length of half long axis of ellips :  A
'Then : half short axis of ellipse is :  B = A sin(lat)
'in which lat is latitude of the place

'x,y coordinates of hourpoints :
'x = A Sin(t)
'y = B Cos(t)
'in which t is the hourangle of the sun corrected for Central Meridian.

'x,y coordinates of scale of date
'x = 0
'y = A cos(lat) tan(decl)
'in which decl is the sun's declination of the date

'x,y of focus of ellipse
'x = A cos(lat)
'y = 0
'And
'x = -A cos(lat)
'y = 0

'An analemmatic sundial isn't very suitable between the 2 tropics
'( Latitude between -23.44 ans 23.44 degrees )

'Coordinates of hourpoints and scale of date are writen to file
'ANALEM3.txt
'Open file in text editor with letter type as Courier or Sanserif

Option Explicit ' Force all variables to be declared before they are used. No adhoc variables

Dim x,y,z,xx,yy,zz,declination,declmax,hourangle,pi As double
Dim latitude,standard_meridian,local_meridian,longitude_correction As Double
Dim flag,count,count1,count2 As Integer
Dim x1,y1,z1,x2,y2 As Double
Dim action,button As String
Dim decl(7),decl2(12),half_long_axis,half_short_axis,halfday,halfdayhour As Double
Dim cx,cy,xa,ya,xb,yb,hrad,wrad As Double
Dim angle,arrow,radius,Time_zone,focus As Double
Dim ys,ye,ms,me,azs,aze As Double
Dim datetext(12),outtext As String
Dim minval, maxval As Double
Dim inputcorrect As boolean
Dim filename As String
Dim text1,text2,text3 As String

dcSetLineParms dcBLACK, dcSOLID, dcTHIN
dcSetCircleParms dcBLACK, dcSOLID, dcTHIN

'Establish the 5 standard line thicknesses in thousands of an inch.
dcSetDrawingData dcLineThin,   .003
dcSetDrawingData dcLineNormal, .008
dcSetDrawingData dcLineThick,  .012
dcSetDrawingData dcLineHeavy,  .024
dcSetDrawingData dcLineWide,   .048

'Maximize the window, close any existing drawing without saving, and start a new drawing.
dcSetDrawingWindowMode dcMaximizeWin
dcCloseWithoutSaving
dcNew ""


'**************************************

'Start of program
init_constants
Input_constants_of_sundial 
half_short_axis = half_long_axis * Sin(Abs(latitude)*zz)
longitude_correction = standard_meridian - local_meridian
halfday = halfdaylength(Abs(latitude),declmax)'in degrees
halfdayhour = halfday/15 'in hours
show_constants 
drawdial
draw_square
drawtext
Write_file
drawsunrise_sunset_circles
dcviewall
'End of program

'**************************************


'Start of subroutines

Sub init_constants
pi = 4 * Atn(1)
zz = pi/180

declmax = 23.44

decl2(1) = "-23.04"
decl2(2) = "-17.23"
decl2(3) = "-7.38"
decl2(4) = "4.75"
decl2(5) = "15.24"
decl2(6) = "22.13"
decl2(7) = "23.07"
decl2(8) = "17.88"
decl2(9) = "8.08"
decl2(10) = "-3.40"
decl2(11) = "-14.60"
decl2(12) = "-21.88"

datetext(1) = "jan"
datetext(2) = "feb"
datetext(3) = "mar"
datetext(4) = "apr"
datetext(5) = "may"
datetext(6) = "jun"
datetext(7) = "jul"
datetext(8) = "aug"
datetext(9) = "sep"
datetext(10) = "oct"
datetext(11) = "nov"
datetext(12) = "dec"
End Sub


Sub Input_constants_of_sundial
Begin Dialog CONSTANTS_INPUT 13,1, 200, 150, "Input constants of sundial"
  Text 32,0,169,10, "HORIZONTAL ANALEMMATIC SUNDIAL"
  Text 32,10,100,10, "with time zone correction" 
  Text 32,32,69,10, "latitude"
  Text 32,48,69,10, "local meridian"
  Text 32,64,69,10, "central meridian"
  Text 32,96,69,10, "half long axis"
  TextBox 132,32,37,10, .latitude
  TextBox 132,48,37,10, .local_meridian
  TextBox 132,64,37,10, .standard_meridian
  TextBox 132,96,37,10, .half_long_axis
  OKButton 84,128,37,12
End Dialog

'Initialize 
Dim prompt As constants_input
prompt.standard_meridian = -15
prompt.local_meridian = -5.5
prompt.latitude = 51.5
prompt.half_long_axis = 100

repeat_until_inputcorrect: 'label to return if input is not correct
action = Dialog(prompt)    'get the input
If test("latitude",prompt.latitude,-90,90) = false Then 
  GoTo repeat_until_inputcorrect
End If
If test("local meridian",prompt.local_meridian,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("central meridian",prompt.standard_meridian,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("half long axis",prompt.half_long_axis,0,1000) = false Then
  GoTo repeat_until_inputcorrect
End If

'Set program variables with input variables, angles in degrees
standard_meridian = prompt.standard_meridian
local_meridian = prompt.local_meridian
latitude = prompt.latitude
half_long_axis = prompt.half_long_axis 
End Sub


Sub show_constants
Dim outtext1,outtext2,outtext3,outtext4 As String
outtext1 = "half short axis = " + Format(half_short_axis,"#####0.00") 
outtext2 = chr$(13) + chr$(13)
outtext3 = "max half day in hours = " + Format(halfdayhour,"#####0.00")
outtext4 = "output in file ANALEM3.txt"
MsgBox outtext1 & outtext2 & outtext3 & outtext2 & outtext4
End Sub


Sub drawdial
'draw ellipse
cx=0
cy=0
If sgn(latitude) <> 0 Then
If sgn(latitude) = 1 Then
xa = half_long_axis*Sin(halfday*zz)
ya = half_short_axis*Cos(halfday*zz)
xb = -xa
yb = ya
End If
If sgn(latitude) = -1 Then
xb = half_long_axis*Sin(halfday*zz)
yb= 0 -half_short_axis*Cos(halfday*zz)
xa = -xb
ya = yb
End If
hrad = half_short_axis
wrad = half_long_axis
dcCreateCircleEx cx,cy,xa,ya,xb,yb,hrad,wrad,0,0
Else
'for latitude 0 draw a line
x1 = half_long_axis
y1 = 0
x2 = 0 -x1
y2 = 0
dcCreateLine x1,y1,x2,y2
End If

'draw part of long axis
x1 = .95 * half_long_axis
y1 = 0
x2 = 0 -x1
y2 = 0
dcCreateLine x1,y1,x2,y2

'draw hourpoints
radius = 10
For count = 0 To 23.999 Step 1
hourangle = (count - 12) * 15 + longitude_correction
If Abs (hourangle) > 180 Then 
hourangle = hourangle - sgn(hourangle)*360
End If
If  Abs(hourangle) <= halfday Then
x = half_long_axis * Sin(hourangle*zz )
y = half_short_axis * Cos (hourangle*zz)*sgn(latitude) 
radius =  half_long_axis/75
dcCreateCircle x,y,radius
End If
Next count

'draw date line for 1th of month
For count = 1 To 6
x1 = 0 - half_long_axis/40
y1 = half_long_axis*Cos(latitude*zz)*Tan(decl2(count)*zz)
x2 = 0
y2 = y1
dcCreateLine x1,y1,x2,y2
Next count
For count = 7 To 12
x1 = half_long_axis/40
y1 = half_long_axis*Cos(latitude*zz)*Tan(decl2(count)*zz)
x2 = 0
y2 = y1
dcCreateLine x1,y1,x2,y2
Next count

'draw centerline
x1 = 0
y1 = half_long_axis*Cos(latitude*zz)*Tan(declmax*zz)
x2 = 0
y2 = 0 -y1
dcCreateLine x1,y1,x2,y2
End Sub


Sub draw_square
x1 = half_long_axis * 1.2
y1 = x1
x2 = 0 - x1
y2 = 0 - x1
dcCreateBox x1,y1,x2,y2
End Sub


Sub drawtext
'Text for dateline
'Draw text only if latitude between -67 and 67 degrees
If (latitude < 67) And (latitude > -67) Then
dcSetTextParms dcRed, "Times New Roman","Standard",0, 4 * half_long_axis,5,0,0
For count = 1 To 6
x1 = 0 - half_long_axis/8
y1 = half_long_axis*Cos(latitude*zz)*Tan(decl2(count)*zz)
outtext = "1-" + datetext(count)
dcCreateText x1,y1,0,outtext
Next count
For count = 7 To 12
x1 = half_long_axis/8
y1 = half_long_axis*Cos(latitude*zz)*Tan(decl2(count)*zz)
outtext = "1-" + datetext(count)
dcCreateText x1,y1,0,outtext
Next count
End If

'Text for hourpoints
dcSetTextParms dcRed, "Times New Roman","Standard",0, 5 * half_long_axis,5,0,0
For count = 0 To 23.999 Step 1
hourangle = (count - 12) * 15 + longitude_correction
If Abs (hourangle) > 180 Then 
hourangle = hourangle - sgn(hourangle)*360
End If
If  Abs(hourangle) <= halfday Then
x1 = half_long_axis * Sin(hourangle*zz )*1.07
y1 = half_short_axis * Cos (hourangle*zz) *1.07*sgn(latitude)
dcCreateText x1,y1,0,CStr(count)
End If
Next count

'Text for N,E,S,W
dcSetTextParms dcBlue, "Times New Roman","Standard",0, 5 * half_long_axis,5,0,0
x1 = 0 - half_long_axis * 1.15
y1 = 0
dcCreateText x1,y1,0,"W"
x1 = half_long_axis * 1.15
y1 = 0
dcCreateText x1,y1,0,"E"

y1 = half_long_axis * 1.15
x1 = 0
dcCreateText x1,y1,0,"N"
y1 = 0 - half_long_axis * 1.15
x1 = 0
dcCreateText x1,y1,0,"S"

'Text for latitude and longitude
dcSetTextParms dcBlack, "Times New Roman","Standard",0, 6 * half_long_axis,5,0,0
outtext = "LAT  "& latitude
x1 = 0 - half_long_axis*.8
y1 = half_long_axis*1.1
dcCreateText x1,y1,0,outtext
outtext = "LONG  "& local_meridian
x1 = half_long_axis*.8
y1 = half_long_axis*1.1
dcCreateText x1,y1,0,outtext
End Sub

Sub Write_file
filename = "analem3.txt"
Open filename For output As #1
Write #1, "ANALEMMATIC SUNDIAL FOR MEAN TIME"
Write #1, "latitude = " & latitude
Write #1, "local_meridian = " & local_meridian 
Write #1, "central meridian = " & standard_meridian
Write #1, chr$(13) 
Write #1, "half long  axis = " & half_long_axis
Write #1, "half short axis = " & Int(half_short_axis*100)/100
focus = half_long_axis * Cos(latitude*zz)
Write #1, "focus = "& Int(focus*100)/100
Write #1, chr$(13) 
Write #1, "HOURPOINTS"
Write #1, "hour       x           y"
For count = 0 To 23.999 Step 1
hourangle = (count - 12) * 15 + longitude_correction
If Abs (hourangle) > 180 Then 
hourangle = hourangle - sgn(hourangle)*360
End If
If  Abs(hourangle) <= halfday Then
x = half_long_axis * Sin(hourangle*zz )
y = half_short_axis * Cos(hourangle* zz)*sgn(latitude)
If count < 10 Then
text1 = " "&CStr(count)
Else
text1 = CStr(count)
End If
x = Int(x*100)/100
y = Int(y*100)/100
text2 = CStr(x)
text3 = CStr(y)
x =12 - Len(text2)
y =12 -Len(text3)
text2 =String(x," ") & text2
text3 =String(y," ") & text3
outtext = text1 & text2 & text3
Write #1, outtext
End If
Next count

Write #1, chr$(13)
Write #1, "SCALE OF DATE"
Write #1, "month      x           y" 
For count = 1 To 12
x = 0 
y = Int(half_long_axis*Cos(latitude*zz)*Tan(decl2(count)*zz)*100)/100
text1 = datetext(count)
text2 = CStr(x)
text3 = CStr(y)
x =11 - Len(text2)
y =12 -Len(text3)
text2 =String(x," ") & text2
text3 =String(y," ") & text3
outtext = text1 & text2 & text3
Write #1,outtext
Next count
Close #1
End Sub


Sub drawsunrise_sunset_circles
dcAddLayer "day_circles"
dcAddLayer "seasonal_markers"
dcSetCurrentLayer "day_circles"
dcSetCircleParms dcGREEN, dcSOLID, dcTHIN
dcSetLineParms dcBLUE, dcSolid, dcTHIN
If latitude <> 0 Then
  latitude = latitude + 0
  'This is a programming trick.
  'Otherwise some daycircles are drawn wrong
  'if latitude is whole number. 
  For count = 1 To 12
  If sgn(decl2(count)) = 1 Then
    focus = half_long_axis * Cos(latitude*zz)
    y = focus * Tan(decl2(count)*zz)
    x = Sqr(focus*focus + y*y)
    radius= Abs(0.5*x/Sin(decl2(count)*zz))
    cy = 0-focus/Tan(2*decl2(count)*zz)
    cx=0
    halfday = halfdaylength(latitude,decl2(count))
    xa = half_long_axis * Sin(halfday*zz )
    ya = half_short_axis * Cos (halfday*zz)*sgn(latitude) 
    xb = 0-xa
    yb = ya
    hrad = radius
    wrad = radius
    dcCreateCircleEx cx,cy,xa,ya,xb,yb,hrad,wrad,0,0
dcSetCurrentLayer "seasonal_markers"
    dcCreateLine 0,y,xa,ya
    dcCreateLine 0,y,xb,yb
dcSetCurrentLayer "day_circles"
  End If
  If sgn(decl2(count)) = -1 Then
    focus = half_long_axis * Cos(latitude*zz)
    y = focus * Tan(decl2(count)*zz)
    x = Sqr(focus*focus + y*y)
    radius= Abs(0.5*x/Sin(decl2(count)*zz))
    cy = 0-focus/Tan(2*decl2(count) *zz)
    cx=0
    halfday = halfdaylength(latitude,decl2(count))
    xb = half_long_axis * Sin(halfday*zz )
    yb = half_short_axis * Cos (halfday*zz)*sgn(latitude) 
    xa = 0-xb
    ya = yb
    hrad = radius
    wrad = radius
    dcCreateCircleEx cx,cy,xa,ya,xb,yb,hrad,wrad,0,0
dcSetCurrentLayer "seasonal_markers"
    dcCreateLine 0,y,xa,ya
    dcCreateLine 0,y,xb,yb
dcSetCurrentLayer "day_circles"
  End If
  Next count
End If
'draw periodic error correction epicycle showing maximum and 
'minimum marker distances at the equinox (e) and solstice (s)
ys = focus * Tan(declmax*zz)
azs = arccos(Sin(declmax*zz)/Cos(latitude*zz))
ms = ys * Tan(azs)
ye = focus * Tan (0.1*zz)
aze = arccos(Sin(0.1*zz)/Cos(latitude*zz))
me = ye * Tan(aze)
radius = (me-ms)/2
x = me - radius
y = 0
dcSetCurrentLayer "seasonal_markers"
dcSetCircleParms dcRed, dcSolid, dcTHIN
dcCreateCircle x,y,radius
x = 0-x
dcCreateCircle x,y,radius
dcSetCurrentLayer "default"
End Sub


Function arcsin(ByVal x) As double
If Abs(x) > 0.999999999999 Then x = sgn(x)*0.999999999999
  arcsin = Atn(x/Sqr(1-x*x))
End Function


Function arccos(ByVal x) As Double
  arccos = pi/2-arcsin(x)
End Function


Function halfdaylength(lat,decl) As Double
  Const maxval = 89.999999999
  Dim var As Double
   
  If Abs(lat-decl) >= maxval Then
    halfdaylength = 0
  Else
    If Abs(lat+decl) >= maxval Then
      halfdaylength = 180
    Else
      var = -Tan(lat*zz)*Tan(decl*zz)
      halfdaylength = arccos(var)/zz
    End If
  End If  
End Function

Function  test(varname,x,minval,maxval) As boolean
If IsNumeric(x) = false Then
test = false
outtext = varname & " must be numeric"
MsgBox outtext
exit Function
End If
If x < minval Or x > maxval Then
outtext = varname & " must be between " & chr$(13) & minval & "  and  " & maxval
MsgBox outtext 
exit Function
End If
test = true
End Function