'**********************************************************************************
'* EoTgraph.bas, a DeltaCad macro to draw EoT graph for a year from -3000 to 4000 *
'*              The year is assumed as in the Gregorian calendar                  * 
'*      Copyright 2000 Fer de Vries & The North American Sundial Society          * 
'*                            http://sundials.org                                 *
'*      This macro may be circulated and modified as long as this header          *
'*                              remains intact.                                   *
'**********************************************************************************
'date 27 june 2000
'date 10 Sep 2009  'Modified by Valentin Hristov lines are commented 'VZH

'REMARK ******************

'The accuracy for sundials is well enough for a man's lifetime.
'It isn't checked what the accuracy will be over some thousands years however the curves
'look very well as in the book Astronomical Algorithms by Jean Meeus, page 174.
'Try year 1246 to see a symmetrical curve.
'In this macro all years are assumed as in the Gregorian calendar, also the years before 1582.

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

'In this Macro formulas dsin, dcos and dtan and reverse functions use degrees.


'MAIN FUNCTIONS AND SUBS IN THIS MACRO :

'dSin(x),  dCos(x)   and dTan(x)   : x is in degrees 
'dAsin(x), dAcos(x), and dAtan(x)  : returns x in degrees

'CalcYearFormula(InitYear)
' --- calculates formulas for EoT and decl. for given year in YearFormula(13)

'Function CalcJulDay(ByVal xyear,xmonth,xday As integer, ByVal xhour As double) As Double
' --- Calculates Julian day out of year,month,day,hour in Gregorian calendar.
'This function is used by CalcYearFormula(InitYear)

'Function EoT(ByVal daynr As double) As double
' --- Calculates Eot out of daynuber
'With variable fact in the procedure output may be choosen to be in
'seconds of time, minutes of time, degrees or radians.

'Function decl(ByVal daynr As double) As Double
' --- Calculates sun's declination out of daynumber
'With variable fact in the procedure output may be choosen to be in
'degrees or radians.

'Function daynr(ByVal y,m,d,h,mi As double) As Double
' --- Calculates daynumber out of year, month, day, hour minutes

'Function halfdayarc(ByVal lat,decl As double) As double
' --- Calculates half day arc out of latitude and declination
'With variable fact in the procedure output may be choosen to be in
'degrees or radians.


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

Dim pi, d2r, r2d As Double
pi = 4*Atn(1)  'Calculates pi
d2r = pi/180   'Factor to convert degrees to radians
r2d = 180/pi   'Factor to convert radians to degrees


'******************************************
'START OF PROGRAM

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

'If calculations of EoT and/or Declination of Sun are needed always start with 
'choosing or inputing InitYear and calculate YearFormula.

Dim InitYear As Integer              '***********
Dim YearFormula(13) As Double        '***********
Dim Longitude,CentralMeridian,LocalCorrDeg As Double  'VZH
Dim LocalCorrMin,CentralLine,EoTTransl,prcorr As Double  'VZH
Dim Sign As Integer

'Input routine
Dim action As String
Dim Normal As Integer
Dim Inverted As Integer

Begin Dialog Year_INPUT 29,58, 131, 100, "Sundial Correction"
  Text 22,2,120,10, "EoT + Longitude Correction"
  Text 32,12,30,10, "Year"
  TextBox 70,12,22,10, .year
  OptionGroup .GROUP_1
    OptionButton 40,22,37,12, "Inverted"
    OptionButton 40,32,41,12, "Normal"
  Text 12,44,40,10, "Longitude"
  TextBox 70,44,40,10, .longitude
  Text 12,56,56,10, "Central meridian"
  TextBox 70,56,40,10, .CentralMeridian
  Text 12,68,115,10, "Decimal deg: East > 0, West < 0"

  OKButton 44,85,37,12
End Dialog

'Initialize 
Dim prompt As Year_input
prompt.year = 2010
prompt.longitude = 23.5011
'prompt.longitude = 35.2261
prompt.CentralMeridian = 30
prompt.GROUP_1 = 1
repeat_until_inputcorrect: 'label to return if input is not correct
action = Dialog(prompt)    'get the input
If test("year",prompt.year,-3000,4000) = false Then 
  GoTo repeat_until_inputcorrect
End If
If test("longitude",prompt.longitude,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("Central Meridian",prompt.CentralMeridian,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
'Set program variables With Input variables

Longitude=prompt.longitude
CentralMeridian=prompt.CentralMeridian

LocalCorrDeg=CentralMeridian-Longitude   'deg
LocalCorrMin=LocalCorrDeg*4              'min
CentralLine=Int((LocalCorrMin+1)/2)*2    'min
EoTTransl=LocalCorrMin-CentralLine       'min

If prompt.GROUP_1 = 0 Then
Sign = 1
End If
If prompt.GROUP_1 = 1 Then
Sign = -1
End If

InitYear = prompt.year               '***********
CalcYearFormula(InitYear)            '***********




Dim EotXY(732) As Double
Dim index,count1,count2,x1,y1,x2,y2,textfact As Double
Dim monthlength(12) As Integer
Dim xmonth,xday,daynum As Double
Dim datetext(12) As String
Dim outtext As String


monthlength(1) = 31
monthlength(2) = 28
If inityear/4 = Int(inityear/4) Then monthlength(2) = 29
monthlength(3) = 31
monthlength(4) = 30
monthlength(5) = 31
monthlength(6) = 30
monthlength(7) = 31
monthlength(8) = 31
monthlength(9) = 30
monthlength(10) = 31
monthlength(11) = 30
monthlength(12) = 31

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"

'draw the grid
textfact = 7
dcSetTextParms dcRed, "Times New Roman","Standard",0, textfact ,5,0,0
'For count1 = -20 To 20 Step 2

For count1 = -20 To 20 Step 2
 prcorr=count1+CentralLine
 If prcorr=int(prcorr/10)*10 Then
  dcSetLineParms dcBLACK, dcSOLID, dcNORMAL
  if prcorr<=0 then
   dcCreateText -53,-Sign*count1*2,0,prcorr
  else
   dcCreateText -53,-Sign*count1*2,0,"+"+CStr(prcorr)
  end if
  if prcorr+60<=0 then
   dcCreateText 53,-Sign*count1*2,0,prcorr+60
  else
   dcCreateText 53,-Sign*count1*2,0,"+"+CStr(prcorr+60)
  end if
 Else
  dcSetLineParms dcBLACK, dcSOLID, dcTHIN
 End If
 dcCreateLine -50,-Sign*count1*2,50,-Sign*count1*2
Next count1

dcSetTextParms dcRed, "Times New Roman","Bold",90, textfact ,21,0,0
dcCreateText -57,0,0,"Civil (WINTER) Time"
dcSetTextParms dcRed, "Times New Roman","Bold",-90, textfact ,21,0,0
dcCreateText 57,0,0,"Daylight Savings (SUMMER) Time"
dcSetTextParms dcRed, "Times New Roman","Bold",0, textfact ,21,0,0
if Longitude=0 then
 dcCreateText 0,-45,0,"Longitude  "+CStr(Longitude)
else
 if Longitude<0 then
  dcCreateText 0,-45,0,"Longitude  "+CStr(abs(Longitude))+" W"
 else
  dcCreateText 0,-45,0,"Longitude  "+CStr(Longitude)+" E"
 end if
end if


dcSetTextParms dcRed, "Times New Roman","Standard",0, textfact ,5,0,0
dcSetLineParms dcBLACK, dcSOLID, dcTHIN

For count1 = 1 To 13
xmonth = count1
daynum = daynr(inityear,xmonth,1,12,0)
x1 = (daynum/2 - 92)/2
y1 = 20*2
x2 = x1
y2 = -y1
dcCreateLine x1,y1,x2,y2
If xmonth < 13 Then
 dcCreateText x1 + 4, -37.7, 0, datetext(xmonth)
 dcCreateText x1 + 4, 37.7, 0, datetext(xmonth)
End If
Next count1

For count1 = -1 To 1 Step 2
x1 = 50*count1
y1 = 20*2
x2 = x1
y2 = -y1
dcCreateLine x1,y1,x2,y2
Next count1

textfact = 12
dcSetTextParms dcRed, "Times New Roman","Standard",0, textfact ,5,0,0
'dcCreateText 0, 38, 0, inityear

textfact = 12
dcSetTextParms dcRed, "Times New Roman","Standard",0, textfact ,5,0,0
'dcCreateText 0, 45, 0, "Equation of Time"
dcCreateText 0, 45, 0, "Correction of the sundial time (in minutes)"



'draw the EoT curve
dcSetLineParms dcBLUE, dcSOLID, dcTHIN
index = 0
For count1 = 1 To 12
xmonth = count1
For count2 = 1 To monthlength(count1) Step 3
xday = count2
daynum = daynr(inityear,xmonth,xday,12,0)
If xmonth = 12 And xday = 31 Then daynum = daynum + 1 'This is 1 january of the next year
'  EotXY(Index+1) = (daynum/2 - 92)/2
'  EotXY(Index+2) = eot(daynum)*2
  EotXY(Index+1) = (daynum/2 - 92)/2
  EotXY(Index+2) = eot(daynum)*2 +EoTTransl*2
index = index + 2
Next count2
dcCreateLine 
Next count1
dcCreateSpline EotXY(1),index/2 ,false   	'draw the EOT


dcCreateBox -60,-50,60,50

dcViewAll				'Make it fit on screen	

'END OF PROGRAM
'******************************************



'******************************************
'The Functions and Sub's

' --- dTan operates on degrees
Function dTan(ByVal value As double) As Double
dTan = Tan(value*d2r)
If Abs(dTan) < 1e-12 Then dTan = 0
End Function

' --- dAtan returns degrees
Function dAtan(ByVal value As double) As Double
dAtan = Atn(value)*r2d
End Function

' --- dSin operates on degrees
Function dSin(ByVal value As double) As Double
dSin = Sin(value * d2r)
If Abs(dSin) < 1e-12 Then dSin = 0
End Function

' --- dAsin returns degrees
Function dAsin(ByVal value As Double) As Double
If Abs(value) > 0.999999999999 Then
dAsin = 90 * sgn(value)
Else
dAsin = dAtan(value/Sqr(1-value*value))
End If
End Function

' --- dCos operates on degrees
Function dCos(ByVal value As double) As Double
dCos = Cos(value * d2r)
If Abs(dCos) < 1e-12 Then dCos = 0
End Function

' --- dAcos returns degrees
Function dAcos(ByVal value As Double) As Double
dAcos = 90 - dAsin(value)
End Function



' --- Calculates Julian day out of year,month,day,hour in Gregorian calendar.
Function CalcJulDay(ByVal xyear,xmonth,xday As integer, ByVal xhour As double) As Double

  Dim help1,help2 As Double
  
  If (xmonth=1) Or (xmonth=2) Then
    xmonth=xmonth+12
    xyear=xyear-1
  End If
  help1=Int(xyear/100)
  help2=2-help1+Int(help1/4)
  CalcJulDay=Int(365.25*xyear)+Int(30.6001*(xmonth+1))+xday+xhour/24+1720994.5+help2
End Function


' --- Calculates a formula for a year to calculate EoT and decl out of daynumber
Sub CalcYearFormula(xyear As Integer)

  Dim l,w,e,epsilon,d,y,fact  As double
  
  fact = r2d*4*60 'convert radians into seconds of time
  'Yearformula 1 - 7   : terms for EoT in seconds of time
  'Yearformula 8       : longitude of sun at 1 jan. 0h:0m:0s
  'Yearformula 9       : obliquity epsilon in degrees
  'Yearformula 10 - 13 : terms for decl. in degrees
  'Calculation of this formula is based on epoch 1900 but for Sundials
  'still valuable. 
  'Literature : 
  'Bulletin of De Zonnewijzerkring, nr. 8, march 1981 and nr. 22, february. 1985,
  'articles by Thijs J. de Vries.

  d = CalcJulDay(xyear,1,0,0)-CalcJulDay(1900,1,0,12)
  l = 279.696678+0.9856473354*d+0.00002267*d*d/100000000.0
  l = l-Int(l/360)*360-360.0
  w = 281.220844+0.0000470684*d+0.00003390*d*d/100000000.0
  w = w-Int(w/360)*360
  e = 0.01675104-0.000011444*d/10000-0.0000000094*d*d/100000000.0
  epsilon= 23.452294-0.0035626*d/10000
  y = dTan(epsilon/2)*dTan(epsilon/2)
  yearformula(1) = (-2*e*dCos(w)-2*e*y*dCos(w))*fact
  yearformula(2) = ( 2*e*dSin(w)-2*e*y*dSin(w))*fact
  yearformula(3) = (y-5/4*e*e*dCos(2*w))*fact
  yearformula(4) = (  5/4*e*e*dSin(2*w))*fact
  yearformula(5) = ( 2*e*y*dCos(w))*fact
  yearformula(6) = (-2*e*y*dSin(w))*fact
  yearformula(7) = (-0.5*y*y)*fact
  yearformula(8) = l
  yearformula(9) = epsilon
  yearformula(10)=  2*e*dCos(w)*r2d
  yearformula(11)= -2*e*dSin(w)*r2d
  yearformula(12)=  5/4*e*e*dCos(2*w)*r2d
  yearformula(13)= -5/4*e*e*dSin(2*w)*r2d
End Sub


' --- Calculates EoT out of daynuber
Function EoT(ByVal daynr As double) As double
'Dim Sign As Integer


  Dim l,help,fact As Double
  'fact = 1          'use this fact for output in seconds of time
  fact = 1/60       'use this fact for output in minutes of time
  'fact = 1/60/4     'use this fact for output in degrees
  'fact = 1/60/4*d2r 'use this fact for output in radians
 
  l=daynr*360/365.2422+yearformula(8)
  help=yearformula(1)*dSin(l)+yearformula(2)*dCos(l)+yearformula(3)*dSin(2*l)
  help=help+yearformula(4)*dCos(2*l)+yearformula(5)*dSin(3*l)
  EoT=help+yearformula(6)*dCos(3*l)+yearformula(7)*dSin(4*l)
  Eot = EoT * fact * Sign
End Function


' --- Calculates sun's declination out of daynumber
Function decl(ByVal daynr As double) As Double

  Dim l,lambda,fact As double
  fact = 1    'use this fact for output in degrees
  'fact = d2r  'use this fact for output in radians
  
  l=daynr*360/365.2422+yearformula(8)
  lambda=l+yearformula(10)*dSin(l)+yearformula(11)*dCos(l)
  lambda=lambda+yearformula(12)*dSin(2*l)+yearformula(13)*dCos(2*l)
  decl=dAsin(dSin(lambda)*dsin(yearformula(9)))
  decl = decl * fact  
End Function


' --- Calculates daynumber out of year, month, day, hour, minutes
Function daynr(ByVal y,M,d,h,mi As double) As Double

  Dim a,help As double

  a=Int((M+9.0)/12.0)
  help=Int(275.0*M/9.0)-2.0*a+d-30.0+h/24.0+mi/60.0/24.0
  If y/4 = Int(y/4) Then
    daynr=help+a
  Else
    daynr=help
  End If
End Function

' --- Calculates half day arc out of latitude and declination
Function halfdayarc(ByVal lat,decl As double) As double

  Const maxval = 89.9999
  Dim help,fact As Double
  fact = 1    'use this fact for output in degrees
  'fact = d2r  'use this fact for output in radians
  
If Abs(lat-decl) >= maxval Then
    halfdayarc = 0
  Else
    If Abs(lat+decl) >= maxval Then
      halfdayarc = 180
    Else
      help = -dTan(lat)*dTan(decl)
      halfdayarc = dAcos(help) * fact
    End If
  End If 
End Function

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