'* SPIDER.BAS
'* Version from February 2008
'* Author: Valentin Hristov
'* E-mail: valhrist@bas.bg
'* Web page: http://www.math.bas.bg/complan/valhrist/mystuff.htm

'* This DeltaCad macro combines the feachures of separate previous
'* files for arbitrary dial plane and also extends them with the
'* new possibility to choose arbitrary directions for the gnomon
'* from the centre.

'* The important particular cases are given as options and the
'* macro calculates the proper direction of the gnomon.

'* Only the periods when the dial plane is illuminated are indicated.
'* If the plane is horisontal, you can find the times of sunrise/sunset.

'* You can switch ON or OFF the Longitude and/or the EOT correction.

'* You can choose how big to be the set of concentric month circles.
'* Giving the radii for 1 Jan and 31 Dec you can make either of them to
'* be inside or outside.

'* In the middle the local time is presented together with the main
'* system of concentric circles for the civil and daylight savings time.

'* There are different layers and you can choose which of them to
'* be ON or OFF for viewing and printing from the "View Layer" menu.
'* In particular, you can keep only one of the two sets of month names,
'* or to remove the local time, etc.

'* Data for the gnomon is given numerically and graphically.

'* There are many important particular cases and the one which gives
'* the name of this macro is the classical spider dial with horizontal
'* dial plane and vertical gnomon for measuring the azimuth of the sun.



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

Dim lat,lon,cm,d,i,rot As Double
Dim pi,d2r,r2d,dm,dmr As Double
Dim w,gnvx,gnvy,gnvz,lgnp,gnpx,gnpy,gnangr,gnang As Double
Dim eot_c,lon_c,feot,flon As Double
Dim x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4 As Double
Dim xx,yy,xxx,yyy,lv,xe,ye,xb,yb,mu,xs,ys,zs,lsg,lsp,gnmax,prmax As Double
Dim lc,lcr,slc,clc,lr,sl,cl,dr,sd,cd,ir,si,ci,rr,sr,cr,lg,lh As Double
Dim rb,re,rmax,rmin,cs,ha,ham,ham1,hac,hacr,shac,chac,dsr,sds,cds As Double

Dim decl(366),eot(366),bm(13),spl(732) As Double

Dim count,button,gnomon_type,Time_int,nl As Double

Dim hhb,hhe,hhv As Boolean

Dim action,p,outtext,datetext(13),lal,lol As String

dcSetLineParms dcBLACK, dcSOLID, dcTHIN
dcSetCircleParms dcBLACK, dcSOLID, dcTHIN
dcSetTextParms dcPurple,"Tahoma","Bold",0,7,21,0,0

'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

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

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

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

Call Sundial

Sub Sundial
'Start of program
init_constants
Input_constants_of_sundial 
If button=0 Then GoTo cncl
Angles
Main
if gnomon_type=0 then LocalTime
dcViewAll
cncl:
End Sub
'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

bm( 1)=  1 '1jan
bm( 2)= 32 '1feb
bm( 3)= 60
bm( 4)= 91
bm( 5)=121
bm( 6)=152
bm( 7)=182
bm( 8)=213
bm( 9)=244
bm(10)=274
bm(11)=305
bm(12)=335 '1dec
bm(13)=365 '31dec
'bm(13)=366 '1jan

For count=1 To 365
w=.017202792*(count-(cm-15)/360)
decl(count)=.367402-23.275*Cos(w+.178044)-.38506*Cos(2*w+.0687076)-.16046*Cos(3*w+.451301)+.00315469*Cos(4*w+.876643)
Next count
decl(366)=decl(1)

For count=1 To 365
w=.017202792*(count-(cm-15)/360)
eot(count)=.00884207-7.36034*Cos(w+1.49487)+9.91068*Cos(2*w-1.21808)-.306344*Cos(3*w+1.82055)+.204774*Cos(4*w-.875115)
Next count
eot(366)=eot(1)

'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"
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,250,196, "Input data for the sundial"

 Text 15,0,300,10, "SPIDER DECLINED AND INCLINED SUNDIAL"
 Text 15,8,150,10, "with corrections for the longitude and the EOT"
 Text 15,20,150,10, "Place"
 TextBox 45,20,122,10, .p
 Text 15,30,150,10, "Latitude (N is >0)"
 TextBox 130,30,37,10, .lat
 Text 15,40,150,10, "Longitude (E is >0)"
 TextBox 130,40,37,10, .lon
 Text 15,50,150,10, "Central meridian"
 TextBox 130,50,37,10, .cm
 Text 15,60,150,10, "Declination of the plane"
 TextBox 130,60,37,10, .d
 Text 15,70,150,10, "  (S=0, W=90, E=-90, N=180 or -180)"
 Text 15,80,150,10, "Inclination of the plane"
 TextBox 130,80,37,10, .i
 Text 15,90,180,10, "  (90 = vertical, from you < 90, towards you >90)"
 Text 15,100,150,10, "Rotation in the plane"
 TextBox 130,100,37,10, .rot
 Text 15,110,180,10, "  (clockwise < 0, anticlockwise >0)"
 Text 15,120,150,10, "Radius for 1 Jan"
 TextBox 130,120,37,10, .rb
 Text 15,130,150,10, "Radius for 31 Dec"
 TextBox 130,130,37,10, .re

 Text 15,140,80,10, "Time interval"
 OptionGroup .time_int
 OptionButton 15,150,35,10, "15 min"
 OptionButton 15,160,35,10, "30 min"
 OptionButton 15,170,35,10, "60 min"

 Text 85,140,80,10, "EOT correction"
 OptionGroup .eot_c
 OptionButton 85,150,30,10, "Yes"
 OptionButton 85,160,30,10, "No"

 Text 155,140,80,10, "Longitude correction"
 OptionGroup .lon_c
 OptionButton 155,150,30,10, "Yes"
 OptionButton 155,160,30,10, "No"

 Text 180,20,80,10, "Gnomon type"
 OptionGroup .gnomon_type
 OptionButton 180,30,80,10, "Polar"
 OptionButton 180,40,80,10, "Vertical"
 OptionButton 180,50,80,10, "Perpendicular"
 OptionButton 180,60,80,10, "Arbitrary with"

 Text 190,70,80,10, "gnomon vector"
 Text 190,80,30,10, "X"
 Text 190,90,30,10, "Y"
 Text 190,100,30,10, "Z"
 TextBox 200,80,37,10, .gnvx
 TextBox 200,90,37,10, .gnvy
 TextBox 200,100,37,10, .gnvz

 OKButton 70,180,37,12
 CANCELButton 140,180,37,12

End Dialog

'Initialize 
Dim prompt As constants_input

prompt.p = "Lozen-Sofia-Bulgaria"
prompt.lat = 42.6
prompt.lon = 23.5
prompt.cm = 30
prompt.d = 180
prompt.i = 0
prompt.rot = 0
prompt.time_int=0
prompt.eot_c=0
prompt.lon_c=0
prompt.rb=1
prompt.re=2.5
prompt.gnomon_type = 0
prompt.gnvx = 1
prompt.gnvy = -1
prompt.gnvz = 1

repeat_until_inputcorrect: 'label to return if input is not correct
button = Dialog(prompt)    'get the input

If (button=0) Then GoTo End_input

If test("Latitude",prompt.lat,-90,90) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("Longitude",prompt.lon,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("Central meridian",prompt.cm,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("Declination",prompt.d,-180,180) = false Then
  GoTo repeat_until_inputcorrect
End If
If test("Inclination",prompt.i,0,180) = false Then
  GoTo repeat_until_inputcorrect
End If

'Set program variables with input variables, angles in degrees
p = prompt.p
lat = CDbl(prompt.lat)
lon = CDbl(prompt.lon)
cm = CDbl(prompt.cm)
d = CDbl(prompt.d)
i = CDbl(prompt.i)
rot = CDbl(prompt.rot)
rb = CDbl(prompt.rb)
re = CDbl(prompt.re)

Time_int=CDbl(prompt.time_int)
If Time_int=0 Then nl=4
If Time_int=1 Then nl=2
If Time_int=2 Then nl=1

eot_c=CDbl(prompt.eot_c)
If eot_c=0 Then feot=1 Else feot=0

lon_c=CDbl(prompt.lon_c)
If lon_c=0 Then flon=1 Else flon=0

gnomon_type=CDbl(prompt.gnomon_type)
If gnomon_type=2 Then
 gnvx=0
 gnvy=0
 gnvz=1
End If
If gnomon_type=3 Then
 gnvx=CDbl(prompt.gnvx)
 gnvy=CDbl(prompt.gnvy)
 gnvz=CDbl(prompt.gnvz)
 If Abs(gnvx)+Abs(gnvy)+Abs(gnvz)=0 Or gnvz<0 Then GoTo repeat_until_inputcorrect
 lg=Sqr(gnvx*gnvx+gnvy*gnvy+gnvz*gnvz)  'length of gnomon
 gnvx=gnvx/lg
 gnvy=gnvy/lg
 gnvz=gnvz/lg                           'length = 1
End If

End_input:

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=CDbl(lat)*d2r   'latitude 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)

rr=rot*d2r   'rotation in radians
sr=Sin(rr)   'sin(rotation)
cr=Cos(rr)   'cos(rotation)

End Sub

Sub Convert

'rotation to horizontal plane
x1=x0
y1=sl*y0+cl*z0
'z1=-cl*y0+sl*z0
'z1=(-1)*cl*y0+sl*z0               '"-cl" in the beginning gives mistake
z1=sl*z0-cl*y0

'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
'z3=(-1)*si*y2+ci*z2
z3=ci*z2-si*y2

'rotation in the plane
x4=cr*x3+sr*y3
'y4=-sr*x3+cr*y3
'y4=(-1)*sr*x3+cr*y3
y4=cr*y3-sr*x3
z4=z3

End Sub

Sub Main

dcAddLayer "Info"
dcAddLayer "Gnomon"
dcAddLayer "Month_Circles"
dcAddLayer "Month_Labels_Up"
dcAddLayer "Month_Labels_Down"
dcAddLayer "Shadow_Casting_Points"
dcAddLayer "Local_Time"

gnmax=0
prmax=0

rmax=(rb+re+Abs(re-rb))/2    ' max(rb,re)
rmin=(rb+re-Abs(re-rb))/2    ' min(rb,re)
dcCreateline -rmax,0,rmax,0
dcCreateline 0,-rmax,0,rmax

If gnomon_type=0 Then
 x0=0           'Celestial N
 y0=0
 z0=1
 convert
 gnvx=x4
 gnvy=y4
 gnvz=z4
End If
If gnomon_type=1 Then
 x0=0           'Vertical
 y0=-cl
 z0=sl
 convert
 gnvx=x4
 gnvy=y4
 gnvz=z4
End If

gnvx=Int(10000000000*gnvx+.5)/10000000000
gnvy=Int(10000000000*gnvy+.5)/10000000000
gnvz=Int(10000000000*gnvz+.5)/10000000000

If gnvz<0 Then      'Celestial S if celestial N is under the plane
 gnvx=-gnvx
 gnvy=-gnvy
 gnvz=-gnvz
End If

For count=1 To 13
 dcSetCurrentLayer "Month_Circles"
 dcCreateCircle 0,0,rb+(re-rb)*(bm(count)-1)/364
 dcSetCurrentLayer "default"
Next count

If Abs(gnvx)+Abs(gnvy)>0 Then
 lgnp=Sqr(gnvx*gnvx+gnvy*gnvy)
 gnpx=gnvx/lgnp
 gnpy=gnvy/lgnp
 dcSetLineParms dcGREEN, dcARROW, dcTHIN
 dcCreateLine rmax*1.1*gnpx,rmax*1.1*gnpy,0,0
 gnangr=Atn(gnvz/lgnp)
 gnang=gnangr*r2d
 dcSetCurrentLayer "Gnomon"
 dcCreateText 0,(0-rmax)*1.5,0,"The angle between the plane "+ _
  "and the gnomon is "+CStr(Int(gnang*10000+.5)/10000)+" degrees."
 dcCreateText 0,(0-rmax)*1.55,0,"The projection of the gnomon "+ _
  "is indicated by an arrow."
 dcCreateLine rmax*.2,(0-rmax)*1.45,0,(0-rmax)*1.45
 dcCreateLine rmax*.2*Cos(gnangr), _
  (0-rmax)*1.45+rmax*.2*Sin(gnangr),0,(0-rmax)*1.45
 dcSetLineParms dcBLACK, dcSOLID, dcTHIN
 dcSetCurrentLayer "default"
Else
 dcSetCurrentLayer "Gnomon"
 dcCreateText 0,(0-rmax)*1.5,0,"The angle between the gnomon "+ _
  "and the plane is "+CStr(90)+" degrees."
 dcCreateLine 0,(0-rmax)*1.45,rmax*.2,(0-rmax)*1.45
 dcCreateLine 0,(0-rmax)*1.45,0,(0-rmax)*1.45+rmax*.2
 dcSetCurrentLayer "default"
End If

dcSetCurrentLayer "Info"
dcSetTextParms dcDarkPurple, "Tahoma","Bold",0,8,21,0,0
dcCreateText 0,rmax*1.18,0,p
If lat>=0 Then lal=CStr(lat)+" N" Else lal=CStr(-Val(lat))+" S"
dcCreateText (0-rmax)*1.2/Sqr(2),rmax*(1.2/Sqr(2)+.1),0,"Latitude "+lal
If lon>=0 Then lol=CStr(lon)+" E" Else lol=CStr(-Val(lon))+" W"
dcCreateText rmax*1.2/Sqr(2),rmax*(1.2/Sqr(2)+.1),0,"Longitude "+lol
dcCreateText (0-rmax)*1.2/Sqr(2),(0-rmax)*(1.2/Sqr(2)+.1),0,"Declination "+CStr(d)
dcCreateText rmax*1.2/Sqr(2),(0-rmax)*(1.2/Sqr(2)+.1),0,"Inclination "+CStr(i)
dcCreateText 0,(0-rmax)*1.18,0,"Rotation "+CStr(rot)
If flon=1 Then
 dcCreateText (0-rmax)*1.3/Sqr(2),(0-rmax)*(1.3/Sqr(2)+.1),0,"With Longitude correction"
Else
 dcCreateText (0-rmax)*1.3/Sqr(2),(0-rmax)*(1.3/Sqr(2)+.1),0,"Without Longitude correction"
End If
If feot=1 Then
 dcCreateText rmax*1.3/Sqr(2),(0-rmax)*(1.3/Sqr(2)+.1),0,"With EOT correction"
Else
 dcCreateText rmax*1.3/Sqr(2),(0-rmax)*(1.3/Sqr(2)+.1),0,"Without EOT correction"
End If
dcSetCurrentLayer "default"

' hour lines
cs=0                         'counter for spline
For ha=0 To 23.999 Step 1/nl
 ' civil and daylight savings labels
 ham=ha-1
 ham=ham-Int(ham/12)*12+1 'ha (12)
 ham1=ham
 ham1=ham1-Int(ham1/12)*12+1 'ha+1 (12)

 hhv=False
 hhb=True
 hhe=False

 For count=1 To 365 Step 2

  hac=ha*15+flon*lc-feot*eot(count)/4
  hacr=hac*d2r
  shac=Sin(hacr)
  chac=Cos(hacr)
  dsr=decl(count)*d2r
  sds=Sin(-dsr)
  cds=Cos(-dsr)
  x0=-shac*cds
  y0=-chac*cds
  z0=sds

  convert

  xx=gnvx-gnvz*x4/z4
  yy=gnvy-gnvz*y4/z4
  lv=Sqr(xx*xx+yy*yy)
  If lv>0 Then
   xx=xx/lv  'unit projection
   yy=yy/lv
   xx=xx*(rb+(re-rb)*(count-1)/364)
   yy=yy*(rb+(re-rb)*(count-1)/364)
  End If

  If z1<0 And z3<0 Then ' illuminated
   cs=cs+1
   spl(2*cs-1)=xx
   spl(2*cs)=yy
   xe=xx
   ye=yy
   If hhb=True Then
    xb=xx
    yb=yy
    hhb=False
    hhv=True
   End If

   If gnvx*z4<>gnvz*x4 Then
    dcSetCurrentLayer "Shadow_Casting_Points"
    mu=z4*xx/(gnvx*z4-gnvz*x4)
    xs=mu*gnvx
    ys=mu*gnvy
    zs=mu*gnvz
    dcSetCircleParms dcGREEN, dcSOLID, dcTHIN
    dcCreateCircle xs,ys,.025
    dcSetLineParms dcBLACK, dcSOLID, dcTHIN
    If Sqr(xs*xs+ys*ys)>prmax Then prmax=Sqr(xs*xs+ys*ys)
    If Sqr(xs*xs+ys*ys+zs*zs)>gnmax Then gnmax=Sqr(xs*xs+ys*ys+zs*zs)
    dcSetCurrentLayer "default"
   End If
  Else
   If cs>2 Then
    If ha=Int(ha) Then
     dcSetLineParms dcBlack, dcSOLID, dcTHICK
     dcSetSplineParms dcBlack, dcSOLID, dcTHICK
    Else
     dcSetLineParms dcBlack, dcSOLID, dcTHIN
     dcSetSplineParms dcBlack, dcSOLID, dcTHIN
    End If
    dcCreateSpline spl(1),cs,False
   Else
    If cs=2 Then
     dcCreateLine spl(1),spl(2),spl(3),spl(4)
    End If
   End If
   cs=0
  End If
 Next count

 hhe=False

 If cs>2 Then                       'if the hour line finishes at 365
  If ha=Int(ha) Then
   dcSetLineParms dcBlack, dcSOLID, dcTHICK
   dcSetSplineParms dcBlack, dcSOLID, dcTHICK
  Else
   dcSetLineParms dcBlack, dcSOLID, dcTHIN
   dcSetSplineParms dcBlack, dcSOLID, dcTHIN
  End If
  dcCreateSpline spl(1),cs,False
 Else
  If cs=2 Then
   dcCreateLine spl(1),spl(2),spl(3),spl(4)
  End If
 End If

 cs=0
 If hhv Then
  If ha=Int(ha) Then
   If rb>re Then
    xx=xb
    yy=yb
    lh=Sqr(xb*xb+yb*yb)
   Else
    xx=xe
    yy=ye
    lh=Sqr(xe*xe+ye*ye)
   End If
   dcSetTextParms dcBlue, "Tahoma","Bold",0,8,21,0,0
   dcCreateText (1+.05*rmax/lh)*xx,(1+.05*rmax/lh)*yy,0,CStr(ham)
   dcSetTextParms dcRed, "Tahoma","Bold",0,8,21,0,0
   dcCreateText (1+.1*rmax/lh)*xx,(1+.1*rmax/lh)*yy,0,CStr(ham1)
  End If
 End If
Next ha

dcSetCurrentLayer "Gnomon"
dcCreateText 0,(0-rmax)*1.65,0,"The needed length ot the gnomon is at least "+ _
  CStr(Int(gnmax*1000000)/1000000)
dcCreateText 0,(0-rmax)*1.7,0,"The corresponding length ot the projection is "+ _
  CStr(Int(prmax*1000000)/1000000)
dcSetCurrentLayer "default"

dcSetTextParms dcPurple,"Tahoma","Bold",0,7,21,0,0
For count=1 To 12
 dcSetCurrentLayer "Month_Labels_Up"
 dcCreateText 0,rb+(re-rb)*(bm(count)+15)/364,0,datetext(count)
 dcSetCurrentLayer "Month_Labels_Down"
 dcCreateText 0,0-(rb+(re-rb)*(bm(count)+15)/364),0,datetext(count)
 dcSetCurrentLayer "default"
Next count

End Sub 'main

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

Sub LocalTime

dcSetCurrentLayer "Local_Time"

dcSetLineParms dcPurple,dcSolid,dcThin

' hour lines
For ha=0 To 23.999 Step 1/nl

 hhb=True

 ham=ha-1
 ham=ham-Int(ham/12)*12+1 'ha (1-12)

 For count=1 To 365 Step 2
  hac=ha*15
  hacr=hac*d2r
  shac=Sin(hacr)
  chac=Cos(hacr)
  dsr=decl(count)*d2r
  sds=Sin(-dsr)
  cds=Cos(-dsr)
  x0=-shac*cds
  y0=-chac*cds
  z0=sds

  convert

  If z1<0 And z3<0 And hhb Then
   hhb=False
   xx=gnvx-gnvz*x4/z4
   yy=gnvy-gnvz*y4/z4
   lv=Sqr(xx*xx+yy*yy)
   If lv>0 Then
    xx=xx/lv  'unit projection
    yy=yy/lv
   End If
   xx=xx*rmin*.8
   yy=yy*rmin*.8
   xxx=xx*.9
   yyy=yy*.9
   If 2*ha=Int(2*ha) Then
    xxx=xx*.8
    yyy=yy*.8
   End If
   If ha=Int(ha) Then
    xxx=xx*.6
    yyy=yy*.6
    dcCreateText xx*1.1,yy*1.1,0,CStr(ham)
   End If
   dcCreateLine xx,yy,xxx,yyy
   GoTo nextha
  End If
 Next count
nextha:
Next ha
dcCreateText 0,.1,0,"LOCAL"
dcCreateText 0,-.1,0,"TIME"
dcSetCurrentLayer "default"

End Sub 'LocalTime

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