'********************************************************
'* NOT FINISHED YET - VARIANT FOR TESTING               *
'* FL.bas is a DeltaCad macro for producing a           *
'* Foster-Lambert Declining Inclining Sundial           *
'* with Longitude Correction and EOT correction         *
'* created by Valentin Hristov (valhrist@bas.bg)        *
'* who was inspired by Mac Oglesby to use               *
'* The North American Sundial Society DeltaCad programs *
'* as tutorials (http://sundials.org)                   *
'********************************************************

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

Dim x,y,z,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3 As Double
Dim l,d,i,lr,dr,ir,sl,cl,sd,cd,si,ci,lon,cm,lc,lcr,slc,clc As Double
Dim xN0,yN0,zN0,xNN0,yNN0,zNN0,lNN0 As Double
Dim xNN,yNN,zNN,XZ0,yZ0,zZ0,xZZ0,yZZ0,zZZ0 As Double
Dim xpr,ypr,zpr,lpr,xpr0,ypr0,zpr0,lpr0,xpr1,ypr1,zpr1,lpr1,xprn,yprn,zprn,lprn As Double
Dim pi,d2r,r2d,h,ch,sh,hr,dm,dmr,dmm,t,g,gr,sg,cg As Double
Dim alpha,beta,betar,gamma, gammar,m,m1,mm,radius As Double
Dim decl(12),eot(13) As Double
Dim datetext(13),outtext,p,ll As String
Dim count,sw,cw As Integer


Dim lat As Double
'Dim ss,dd,ll As Integer
'Dim lat,dw,iw,ha,declmax As Double
'Dim alpha,beta,gamma,v,b,ts,M,mm,m1,radius,nN,cnN,snN,sN As Double
'Dim radius,flag,count,count1,count2,ss As Integer
Dim action,button As String
Dim cx,cy,xa,ya,xb,yb,hrad,wrad As Double
Dim inputcorrect As boolean
Dim filename As String
Dim text1,text2,text3 As String
Dim dx,dy As Double

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 
Angles
Vectors
Preparation
DateLine
Hours
EquationOfTime
Control
'End of program

'**************************************
'Start of subroutines

Sub init_constants
pi = 4 * Atn(1)
d2r = pi/180
r2d = 180/pi
dm = 23.43954
dmr = dm*d2r

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

eot(1) = "3.4"
eot(2) = "13.6"
eot(3) = "12.5"
eot(4) = "4.1"
eot(5) = "-2.8"
eot(6) = "-2.3"
eot(7) = "3.6"
eot(8) = "6.3"
eot(9) = "0.2"
eot(10) = "-10.1"
eot(11) = "-16.3"
eot(12) = "-11.2"
eot(13) = "3.4"

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"
datetext(13) = "jan"
End Sub

Sub Input_constants_of_sundial
Begin Dialog CONSTANTS_INPUT 13,1,200,130, "Input data for the sundial"
 Text 15,0,300,10, "FOSTER-LAMBERT DECLINED INCLINED SUNDIAL"
 Text 15,8,150,10, "with corrections for the latitude and the EOT"
 Text 15,20,150,10, "Place"
 Text 15,32,150,10, "Latitude (N is >0)"
 Text 15,44,150,10, "Longitude (E is >0)"
 Text 15,56,150,10, "Central meridian"
 Text 15,68,150,10, "Declination of the plane"
 Text 15,80,150,10, "   (S=0, W=90, E=-90, N=180 or -180)"
 Text 15,92,150,10, "Inclination of the plane"
 Text 15,104,180,10, "   (90 - vertical, from you < 90, towards you >90)"
 TextBox 88,20,99,10, .p
 TextBox 150,32,37,10, .l
 TextBox 150,44,37,10, .lon
 TextBox 150,56,37,10, .cm
 TextBox 150,68,37,10, .d
 TextBox 150,92,37,10, .i
 OKButton 82,116,37,12
End Dialog

'Initialize 
Dim prompt As constants_input

'prompt.p = "Blattleboro - USA"
'prompt.l = 42.85
'prompt.lon = -72.55
'prompt.cm = -75
'prompt.d = -46.56
'prompt.i = 90

prompt.p = "Lozen - Sofia - Bulgaria"
prompt.l = 42.6
prompt.lon = 23.5
prompt.cm = 30
prompt.d = 5.5
prompt.i = 90

'prompt.p = "Test variables"
'prompt.l = 40
'prompt.lon = 26.25
'prompt.cm = 30
'prompt.d = 0
'prompt.i = 90

repeat_until_inputcorrect: 'label to return if input is not correct
action = Dialog(prompt)    'get the input
If test("l",prompt.l,-90,90) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("lon",prompt.lon,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("cm",prompt.cm,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("d",prompt.d,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("i",prompt.i,0,180) = false Then
  GoTo repeat_until_inputcorrect
End If

'Set program variables with input variables, angles in degrees
p = prompt.p
l = prompt.l
lon = prompt.lon
cm = prompt.cm
d = prompt.d
i = prompt.i
End Sub

Sub Angles
lc=lon-cm    'longitude correction
lcr=lc*d2r   'longitude correction in radians
slc=sin(lcr) 'sin(longitude correction)
clc=cos(lcr) 'cos(longitude correction)
lr=l*d2r     'longitude in radians
sl=sin(lr)   'sin(latitude)
cl=cos(lr)   'cos(latitude)
dr=d*d2r     'declination in radians
sd=sin(dr)   'sin(declination)
cd=cos(dr)   'cos(declination)
ir=i*d2r     'inclination in radians
si=sin(ir)   'sin(inclination)
ci=cos(ir)   'cos(inclination)
End Sub

Sub Vectors
'Celestial North
xN0=0
yN0=0
zN0=1
' 0 hours without longitude correction
xZ0=0
yZ0=-1
zZ0=0
' 0 hours with longitude correction
xZZ0=-slc
yZZ0=-clc
zZZ0=0
End Sub

Sub Preparation
dcCreateCircle 0,0,.02

x0=xN0
y0=yN0
z0=zN0
convert
xNN0=x3                       'Celestian N in new coordinates
yNN0=y3
zNN0=z3
lNN0=Sqr(xNN0*xNN0+yNN0*yNN0)
xNN=xNN0/lNN0                 'Unit projection of N in the plane
yNN=yNN0/lNN0
if zNN0<0 then cw=-1 else cw=1 'Clockwise flag
if cw=1 then
 xpr=xNN0
 ypr=yNN0
 zpr=1+zNN0
else
 xpr=-xNN0
 ypr=-yNN0
 zpr=1-zNN0
end if
lpr0=Sqr(xpr*xpr+ypr*ypr+zpr*zpr)
xpr0=xpr/lpr0                    'Unit vector for projection
ypr0=ypr/lpr0
zpr0=zpr/lpr0
lpr1=Sqr(xpr*xpr+ypr*ypr)
xpr1=xpr/lpr1                    'Unit vector in the plane
ypr1=ypr/lpr1
zpr1=0
if xpr1>=0 then
 betar=arcsin(ypr1)
else
 if ypr1>=0 then
  betar=arccos(xpr1)
 else
  betar=-arccos(xpr1)
 end if
end if
beta=betar*r2d  'angle ot the projection in the plane

cg=xpr0*xpr1+ypr0*ypr1+zpr0*zpr1 'cos(gnomon)
gr=arccos(cg) 'angle of gnomon in radians
g=gr*r2d      'angle of gnomon in degrees
sg=Sqr(1-cg*cg) 'sin(gnomon)
End Sub

Sub DateLine
'draw the date line for 1st of month
alpha=180-2*g
M=Tan(alpha/2*d2r)
m1=2*M*Tan(dmr)
dcCreateLine 0,m1,0,-m1
dcCreateLine 0,1.8,0,2.2
dcCreateLine 0,-1.8,0,-2.2
dcCreateLine 1.8,0,2.2,0
dcCreateLine -1.8,0,-2.2,0
dcSetTextParms dcDarkgreen, "Times New Roman","Standard",0,6,4,0,0
'date marks
For count = 1 To 6
x1 = .07
y1 = 2*M*Tan(decl(count)*d2r)
x2 = 0
y2 = y1
dcCreateLine x1,y1,x2,y2
x1 = .1
outtext = "1-" + datetext(count)
dcCreateText x1,y1,0,outtext
Next count
dmm=2*M*Tan(dmr)
dcCreateText -.04,dmm+.06,0,"21 jun"
dcCreateLine -.02,dmm,.02,dmm
dcSetTextParms dcDarkgreen, "Times New Roman","Standard",180,6,4,0,0
For count = 7 To 12
x1 = 0-.07
y1 = 2*M*Tan(decl(count)*d2r)
x2 = 0
y2 = y1
dcCreateLine x1,y1,x2,y2
x1 = 0-.1
outtext = "1-" + datetext(count)
dcCreateText x1,y1,0,outtext
Next count
dcCreateText .04,-dmm-.06,0,"23 dec"
dcCreateLine -.02,-dmm,.02,-dmm
dcSelectObjInBox -2,-2,2,2
dcSetSelectBase 0,0
dcRotateSelObjs (-90*cw+beta)
dcUnSelectAll
End Sub

Sub Hours
dcCreateLine 0,1.8,0,2.2
dcCreateLine 0,-1.8,0,-2.2
dcCreateLine 1.8,0,2.2,0
dcCreateLine -1.8,0,-2.2,0
dcCreateLine -.05,1.9,0,1.8
dcCreateLine 0,1.8,.05,1.9
'dcCreateLine 0,-2,0,2
'dcCreateLine -2,0,2,0
x0=xZZ0
y0=yZZ0
z0=zZZ0
convert
t=-z3/zpr0
x=x3+xpr0*t
y=y3+ypr0*t
if x>=0 then
 gammar=arcsin(y)
else
 if y>=0 then
  gammar=arccos(x)
 else
  gammar=-arccos(x)
 end if
end if
gamma=gammar*r2d

dcCreateCircle 2*x,2*y,.07

dcCreateCircle 0,0,2.2
dcCreateCircle 0,0,2
dcCreateCircle 0,0,1.8

for count=0 to 23.999 step 1
' h=count*15*d2r
' sh=sin(h)
' ch=cos(h)
' x0=ch*xZZ0+sh*yZZ0
' y0=-sh*xZZ0+ch*yZZ0
' z0=0
' convert
' t=-z3/zPRn
' x=x3+xPRn*t
' y=y3+yPRn*t
' dcCreateCircle 2*x,2*y,.02
h=-(count*15)*cw+gamma
hr=h*d2r
sh=sin(hr)
ch=cos(hr)
x=2*ch
y=2*sh
dcCreateCircle x,y,.04
dcSetTextParms dcBlue, "Times New Roman","Standard",90+h,10,5,0,0
dcCreateText x*.94,y*.94,0,CStr(count)
dcSetTextParms dcRed, "Times New Roman","Standard",-90+h,10,5,0,0
dcCreateText x*1.06,y*1.06,0,CStr(count+1)
next count
for count=0 to 23.999 step .5
h=-(count*15)*cw+gamma
hr=h*d2r
sh=sin(hr)
ch=cos(hr)
x=2*ch
y=2*sh
dcCreateCircle x,y,.03
next count
for count=0 to 23.999 step .25
h=-(count*15)*cw+gamma
hr=h*d2r
sh=sin(hr)
ch=cos(hr)
x=2*ch
y=2*sh
dcCreateCircle x,y,.02
next count

'dcCreateCircle 2*xpr1,2*ypr1,.06
'drawing the gnomon
dcCreateLine 0,0,cg*xpr1+sg*ypr1,-sg*xpr1+cg*ypr1
End Sub

Sub EquationOfTime
For M=-15 To 15 Step 5
mm=(90+M/4)*d2r
x=Cos(mm)
y=Sin(mm)
x1=1.75*x
y1=1.75*y
dcCreateLine x,y,x1,y1
Next M
For count=1 To 13
radius=1.78-.06*count
x=radius*Cos((90+cw*eot(count)/4)*d2r)
y=radius*Sin((90+cw*eot(count)/4)*d2r)
dcCreateCircle x,y,.01
dcCreateCircle x,y,.02
dcCreateCircle x,y,.03
dcSetTextParms dcRed, "Times New Roman","Standard",0,6,5,0,0
dcCreateText x+.11,y,0,datetext(count)
dcCreateText x-.07,y,0,"1"
Next fcount
End Sub

Sub Control
dcSetTextParms dcPurple, "Times New Roman","Standard",0,12,21,0,0
dcCreateText 0,-1.25,0,p
if l<0 then ll="S" else ll="N"
dcCreateText -.8,1.2,0,CStr(l)+" "+ll
if lon<0 then ll="W" else ll="E"
dcCreateText .8,1.2,0,CStr(lon)+" "+ll
dcSetTextParms dcRed, "Times New Roman","Standard",0,10,24,0,0
dcCreateText 2.25,2   ,0,"latitude="+CStr(l)
dcCreateText 2.25,1.85,0,"longitude="+CStr(lon)
dcCreateText 2.25,1.7 ,0,"central meridian=" +CStr(cm)
dcCreateText 2.25,1.55,0,"longitude correction=" +CStr(lc)
dcCreateText 2.25,1.4 ,0,"declination=" +CStr(d)
dcCreateText 2.25,1.25,0,"inclination=" +CStr(i)
dcCreateText 2.25,1.1 ,0,"Celestial North=("+CStr(xN0)+","+CStr(yN0)+","+CStr(zN0)+")"
dcCreateText 2.25,0.95,0,"New North=("+CStr(xNN0)+","+CStr(yNN0)+","+CStr(zNN0)+")"
dcCreateText 2.25,0.8 ,0,"0 hours=("+CStr(xZZ0)+","+CStr(yZZ0)+","+CStr(zZZ0)+")"
dcCreateText 2.25,0.65,0,"lNN0="+CStr(lNN0)
dcCreateText 2.25,0.5 ,0,"lpr0="+CStr(lpr0)
dcCreateText 2.25,.35 ,0,"lPRn="+CStr(lPRn)
dcCreateText 2.25,.2  ,0,"g="+CStr(g)
dcCreateText 2.25,.05 ,0,"alpha="+CStr(alpha)
dcCreateText 2.25,-.1 ,0,"beta="+CStr(beta)
dcCreateText 2.25,-.25,0,"gamma="+CStr(gamma)
'dcCreateText 2.25,-.4 ,0,"sN="+CStr(sN)
End Sub

Sub Convert
'rotation to horizontal plane
x1=x0
y1=sl*y0+cl*z0
z1=-cl*y0+sl*z0
'rotation for declination of the plane
x2=cd*x1-sd*y1
y2=sd*x1+cd*y1
z2=z1
'rotation for the inclination of the plane
x3=x2
y3=ci*y2+si*z2
z3=-si*y2+ci*z2
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 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