'********************************************************
'* NOT FINISHED YET - VARIANT FOR TESTING               *
'* SDBIFARB.bas is a DeltaCad macro for producing a     *
'* Bifilar Sundial with two arbitrary gnomons and       *
'* arbitrary dial plane with Longitude Correction and   *
'* EOT Correction                                       *
'* Created by Valentin Hristov (valhrist@bas.bg) in the *
'* end of December 2007 as New Year Present.            *
'* (www.math.bas.bg/complan/valhrist/mystuff.htm)       *
'* I was inspired by Mac Oglesby to use the             *
'* North American Sundial Society DeltaCad programs     *
'* as tutorials (http://sundials.org)                   *
'* and made with DeltaCad (www.deltacad.com)            *
'* different types of sundials.                         *
'********************************************************

'* The classical bifilar sundial has two gnomons which are perpendicular
'* to each other and parallel to the plane with projections from top to
'* bottom and from left to right (as in my previous macro).
'* This DeltaCad macro allows mach biger freedom for choice of:
'*   - dial plane
'*   - two arbitrary straight gnomons given by:
'*      1) initial point
'*      2) directional vector
'* In addition to this it is possible to choose the option for drawing a
'* Double Bifilar Dial with Left-Right symmetry.
'* I saw a picture of such dial on the page
'*   http://www.de-zonnewijzerkring.nl/summaries/home-summaries-2006-3.htm
'* and this encouraged me to add much more flexibility to my previous
'* macro (December 2007) for classical bifilar dial.

'* Now you can test "almost infinitely" many possibilities.
'* I would suggest to change slightly the initial Y-coordinate 2 of any
'* of the gnomons to 1.9 and 2.1 and to see how the hour lines rotate
'* remaining parallel to each other.
'* Another test could be to replace the X-coordinate 0 of the second
'* gnomon to -.4 and to see all curved hour lines except 12 hour.
'* It is also interesting to see the result if the EoT correction and/or
'* the Longitude correction are/is included.

'* Different layers are drawn when you run the macro but I decided to
'* hide some of them in the end. To see them you have to switch them ON
'* in the "View-Layers" menu.

'* Before printing switch OFF the not needed layers.

'* You can draw as particular case the classical bifilar dial by
'* choosing the first gnomon with initial point (0,0,1) and vector
'* (0,1,0). The second gnomon should be with (0,0,h) and (1,0,0) resp.
'* If the height h is suitably choosen, then the hour lines could be
'* equiangular (15 degrees per hour). For Decl=0, Incl=0, Rot=0, the
'* height should be SIN(latitude). You can give symmetric boundaries for
'* the printed area.

'* In the particular case of common point of the gnomons the usual NODUS
'* DIAL is drawn with arbitrary position of the dial plate.
'* Then the gnomons are not needed and only a shadow casting point above
'* the plane at the given height is enough.

'* Some parts of the macro code are not needed but I do not remove them
'* and only make the lines activating them commented (using ').

'* You can see the picture of another type of sundial generated by one
'* of my DeltaCad macros at
'* www.flickr.com/photos/Valentin_Hristov/261303801/
'* Click on the button "All sizes" to see a bigger photo with details.
'* I am very grateful to my friends Daniela (www.danyo.net) and
'* Todor (www.todor.org) who converted my drawing into a real art piece!!!

'* There is a Demo version of DeltaCad at www.deltacad.com.
'* Use the menu "Options - Macro - Run..." or the separate "Macro"
'* button - "Edit Macro List", add the file, and "Run Macro".
'* In the dialog box use DECIMAL DEGREES. Negative values indicate
'* South for latitude and West for longitude and central meridian.

'* E N J O Y !!!

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

dcSetLineParms dcBlack, dcSolid, dcThin

Dim l,p,lon,cm,rot,i,pi,d2r,r2d,dm,dmr,w,ll,ud,lol,s,cs As Double
Dim ha,ham,ham1,has,has1,lc,hac,tir,hacr,thac,xe,ye,xb,yb As Double
Dim x,y,xx,yy,d,h,lcr,slc,clc,lr,sl,cl,dr,sd,cd,ir,si,ci As Double
Dim xN0,yN0,zN0,xZ0,yZ0,zZ0,xZZ0,yZZ0,zZZ0,rr,sr,cr as Double
Dim xNN0,yNN0,zNN0,lNN0,xNN,yNN,xpr,ypr,zpr,dsr,sds,cds as Double
Dim lpr0,xpr0,ypr0,zpr0,lpr1,xpr1,ypr1,zpr1,cg,gr,g,sg as Double
Dim x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,shac,chac as Double
Dim beta,betar,sb,cb,feot,flon,dbl,fdbl,ceot,clon as Double
Dim r,time_int,eot_c,lon_c,per,spring,bb,ee,a,b,at,ab,al,ar,fh as Double
Dim g1px,g1py,g1pz,g1vx,g1vy,g1vz as Double
Dim g2px,g2py,g2pz,g2vx,g2vy,g2vz as Double
Dim g1x,g1y,g1z,g2x,g2y,g2z,l1,l2,m1,m2 as Double
Dim g1nx,g1ny,g1nz,g2nx,g2ny,g2nz as Double

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

Dim button,count,nl,bm(13),ind(37),nom,ind1 As Integer

Dim action,outtext As String

Dim datetext(13) As String

Dim hhv,hhb,hhe,nlcheck As Boolean


dcSetLineParms dcBLACK, dcSOLID, dcTHIN
dcSetCircleParms dcBLACK, dcSOLID, dcTHIN
dcSetTextParms dcPurple, "Tahoma","Bold",0,12,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

'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
Main
'Latitude
'End of program


'**************************************
'Start of subroutines
'''''''''''''''''''''''''''''''''''''''
Sub Input_constants_of_sundial

Begin Dialog CONSTANTS_INPUT 10,25,295,206, "Input data for the sundial"
 Text 15,0,180,10, "BIFILAR SUNDIAL WITH ARBITRARY ORIENTATION"
 Text 15,10,180,10, "with corrections for the latitude and the EOT"
 Text 15,20,150,10, "Place"
 TextBox 88,20,99,10, .p
 Text 15,30,180,10, "Latitude(N>0)  Longitude(E>0)  Central meridian(E>0)"
 TextBox 15,40,30,10, .l
 TextBox 65,40,30,10, .lon
 TextBox 115,40,30,10, .cm
 Text 15,50,150,10, "  Use decimal degrees !!!"
 Text 15,60,150,10, "Declination of the dial plane"
 TextBox 150,60,37,10, .d
 Text 15,70,150,10, "   (S=0, W>0, E<0, N=180 or -180)"
 Text 15,80,150,10, "Inclination of the dial plane"
 TextBox 150,80,37,10, .i
 Text 15,90,180,10, "   (vertical = 90, from you < 90, towards you >90)"
 Text 15,100,150,10, "Rotation in the dial plane"
 TextBox 150,100,37,10, .rot
 Text 15,110,180,10, "   (vertical = 0, clockwise < 0, anticlockwise > 0)"

 Text 15,130,50,10, "Period"
 OptionGroup .per
 OptionButton 38,120,85,10, "21 Dec - 21 Jun"
 OptionButton 38,130,85,10, "21 Jun - 21 Dec"
 OptionButton 38,140,85,10, "Whole year"

 Text 105,130,80,10, "Time interval"
 OptionGroup .time_int
 OptionButton 150,120,35,10, "15 min"
 OptionButton 150,130,35,10, "30 min"
 OptionButton 150,140,35,10, "60 min"

 Text 15,155,80,10, "EOT correction"
 OptionGroup .eot_c
 OptionButton 65,150,30,10, "Yes"
 OptionButton 65,160,30,10, "No"

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

' Text 15,170,180,10, "Heights of the gnomons:   Up<->Down     Left<->Right"
' TextBox 100,180,37,10, .a
' TextBox 150,180,37,10, .b

 Text 15,170,180,10, "Rectangular area: Bottom / Top / Left / Right"
 TextBox 15,180,37,10, .ab
 TextBox 60,180,37,10, .at
 TextBox 105,180,37,10, .al
 TextBox 150,180,37,10, .ar

 Text 200,20,37,10, "Gnomon 1"
 Text 200,30,80,10, "Point            Vector"
 TextBox 200,40,37,10, .g1px
 TextBox 200,50,37,10, .g1py
 TextBox 200,60,37,10, .g1pz
 Text 239,40,37,10, "X"
 Text 239,50,37,10, "Y"
 Text 239,60,37,10, "Z"
 TextBox 246,40,37,10, .g1vx
 TextBox 246,50,37,10, .g1vy
 TextBox 246,60,37,10, .g1vz
 Text 200,80,37,10, "Gnomon 2"
 Text 200,90,80,10, "Point            Vector"
 TextBox 200,100,37,10, .g2px
 TextBox 200,110,37,10, .g2py
 TextBox 200,120,37,10, .g2pz
 Text 239,100,37,10, "X"
 Text 239,110,37,10, "Y"
 Text 239,120,37,10, "Z"
 TextBox 246,100,37,10, .g2vx
 TextBox 246,110,37,10, .g2vy
 TextBox 246,120,37,10, .g2vz

 Text 200,155,80,10, "Double bifilar"
 OptionGroup .dbl
 OptionButton 250,150,30,10, "Yes"
 OptionButton 250,160,30,10, "No"
 Text 200,170,150,10, "(Left<->Right symmetry)"

 OKButton 42,192,37,12
 CANCELButton 116,192,37,12
End Dialog

'Initialize 
Dim prompt As constants_input

prompt.p = "Lozen-Sofia-Bulgaria"
prompt.l = 42.6
prompt.lon = 23.5
prompt.cm = 30
prompt.d = 0
prompt.i = 42.6
prompt.rot = 0
prompt.per = 0
prompt.time_int = 0
prompt.eot_c = 1
prompt.lon_c = 1
'prompt.a = 1
'prompt.b = 0.676876
prompt.ab = -1
prompt.at = 10
prompt.al = -1
prompt.ar = 10
prompt.g1px = 0
prompt.g1py = 0
prompt.g1pz = 0
prompt.g1vx = 0
prompt.g1vy = 2
prompt.g1vz = 1
prompt.g2px = 0
prompt.g2py = 1
prompt.g2pz = 0
prompt.g2vx = .25
prompt.g2vy = 2
prompt.g2vz = 1
prompt.dbl = 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("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
If test("rot",prompt.rot,-90,90) = false Then
  GoTo repeat_until_inputcorrect
End If
'If test("r",prompt.r,0,10) = false Then
'  GoTo repeat_until_inputcorrect
'End If

'Set program variables with input variables, angles in degrees
p = prompt.p
l = CDbl(prompt.l)
lon = CDbl(prompt.lon)
cm = CDbl(prompt.cm)
d = CDbl(prompt.d)
i = CDbl(prompt.i)
rot = CDbl(prompt.rot)
'a = CDbl(prompt.a)
'b = CDbl(prompt.b)
at = CDbl(prompt.at)
ab = CDbl(prompt.ab)
al = CDbl(prompt.al)
ar = CDbl(prompt.ar)
per = prompt.per
time_int=prompt.time_int
eot_c=prompt.eot_c
lon_c=prompt.lon_c
g1px=CDbl(prompt.g1px)
g1py=CDbl(prompt.g1py)
g1pz=CDbl(prompt.g1pz)
g1vx=CDbl(prompt.g1vx)
g1vy=CDbl(prompt.g1vy)
g1vz=CDbl(prompt.g1vz)
g2px=CDbl(prompt.g2px)
g2py=CDbl(prompt.g2py)
g2pz=CDbl(prompt.g2pz)
g2vx=CDbl(prompt.g2vx)
g2vy=CDbl(prompt.g2vy)
g2vz=CDbl(prompt.g2vz)
dbl=prompt.dbl
if time_int=0 then nl=4
if time_int=1 then nl=2
if time_int=2 then nl=1
if eot_c=0 then feot=1 else feot=0
if lon_c=0 then flon=1 else flon=0
if at<=ab or at<0 or ab>0 then GoTo repeat_until_inputcorrect
if ar<=al or al>0 or ar<0 then GoTo repeat_until_inputcorrect
if abs(g1vx)+abs(g1vy)+abs(g1vz)=0 then GoTo repeat_until_inputcorrect
if abs(g2vx)+abs(g2vy)+abs(g2vz)=0 then GoTo repeat_until_inputcorrect
if dbl=0 then fdbl=1 else fdbl=0

end_input:

End Sub

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

ind( 1)=355
ind( 2)=  1 '1Jan
ind( 3)= 11
ind( 4)= 21
ind( 5)= 32 '1Feb
ind( 6)= 42
ind( 7)= 52
ind( 8)= 60 '1Mar
ind( 9)= 70
ind(10)= 80
ind(11)= 91 '1Apr
ind(12)=101
ind(13)=111
ind(14)=121 '1May
ind(15)=131
ind(16)=141
ind(17)=152 '1Jun
ind(18)=162
ind(19)=172
ind(20)=182 '1Jul
ind(21)=192
ind(22)=202
ind(23)=213 '1Aug
ind(24)=223
ind(25)=233
ind(26)=244 '1Sep
ind(27)=254
ind(28)=264
ind(29)=274 '1Oct
ind(30)=284
ind(31)=294
ind(32)=305 '1Nov
ind(33)=315
ind(34)=325
ind(35)=335 '1Dec
ind(36)=345
ind(37)=355

'ind( 1)=  1 '1Jan
'ind( 2)= 11
'ind( 3)= 21
'ind( 4)= 32 '1Feb
'ind( 5)= 42
'ind( 6)= 52
'ind( 7)= 60 '1Mar
'ind( 8)= 70
'ind( 9)= 80
'ind(10)= 91 '1Apr
'ind(11)=101
'ind(12)=111
'ind(13)=121 '1May
'ind(14)=131
'ind(15)=141
'ind(16)=152 '1Jun
'ind(17)=162
'ind(18)=172
'ind(19)=182 '1Jul
'ind(20)=192
'ind(21)=202
'ind(22)=213 '1Aug
'ind(23)=223
'ind(24)=233
'ind(25)=244 '1Sep
'ind(26)=254
'ind(27)=264
'ind(28)=274 '1Oct
'ind(29)=284
'ind(30)=294
'ind(31)=305 '1Nov
'ind(32)=315
'ind(33)=325
'ind(34)=335 '1Dec
'ind(35)=345
'ind(36)=355
'ind(37)=366 '1Jan

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)=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) = "J"
'datetext( 2) = "F"
'datetext( 3) = "M"
'datetext( 4) = "A"
'datetext( 5) = "M"
'datetext( 6) = "J"
'datetext( 7) = "J"
'datetext( 8) = "A"
'datetext( 9) = "S"
'datetext(10) = "O"
'datetext(11) = "N"
'datetext(12) = "D"
'datetext(13) = "J"

'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 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)

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

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=x4                       'Celestian N in new coordinates
yNN0=y4
zNN0=z4
lNN0=Sqr(xNN0*xNN0+yNN0*yNN0)
xNN=xNN0/lNN0                 'Unit projection of N in the plane
yNN=yNN0/lNN0
if zNN0<0 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 of the projection in the plane
sb=sin(betar)
cb=cos(betar)

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)
'dcCreateLine 0,0,cg,sg
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
'rotation of the box in its plane
x4=cr*x3+sr*y3
y4=-sr*x3+cr*y3
z4=z3
End Sub

Sub Main


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
dcAddLayer "ShadowCastingPoints"
dcAddlayer "DateLines"
dcAddlayer "Gnomons"
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
if fdbl=1 then
 if abs(al)>abs(ar) then
  ar=-abs(al)
 else
  al=-abs(ar)
 end if
end if

dcCreateLine 0,ab,0,at
dcCreateLine al,0,ar,0
dcCreateLine al,ab,al,at
dcCreateLine al,at,ar,at
dcCreateLine ar,at,ar,ab
dcCreateLine ar,ab,al,ab

'dcSetCurrentLayer "ShadowCastingPoints"

'dcSetLineParms dcBlack, dcStitch, dcThin

'dcCreateLine 0,ab-a,0,ab
'dcCreateLine 0,ab-a,0-a,ab
'dcCreateLine 0,ab-a,a,ab

'dcCreateLine 0,at+a,0,at
'dcCreateLine 0,at+a,0-a,at
'dcCreateLine 0,at+a,a,at

'dcCreateLine al-b,0,al,0
'dcCreateLine al-b,0,al,0-b
'dcCreateLine al-b,0,al,b

'dcCreateLine ar+b,0,ar,0
'dcCreateLine ar+b,0,ar,0-b
'dcCreateLine ar+b,0,ar,b

'dcSetCurrentLayer "default"

'dcCreateCircle 0,ab-a,.05
'dcCreateCircle 0,at+a,.05
'dcCreateCircle al-b,0,.05
'dcCreateCircle ar+b,0,.05

if l>=0 then
 ll=CStr(l)+" N"
 ud=0
else
 ll=CStr(-Val(l))+" S"
 ud=180
end if
if lon>=0 then lol=CStr(lon)+" E" else lol=CStr(-Val(lon))+" W"
dcSetTextParms dcDarkPurple, "Tahoma","Bold",0,8,21,0,0
dcCreateText (al+ar)/2,at+0.2,0,p
dcCreateText (al+ar)/2-2.5,at+0.2,0,"Latitude "+ll
dcCreateText (al+ar)/2+2.5,at+0.2,0,"Longitude "+lol
dcCreateText (al+ar)/2+1.5,ab-0.2,0,"Declination "+CStr(d)
dcCreateText (al+ar)/2-1.5,ab-0.2,0,"Inclination "+CStr(i)
dcCreateText (al+ar)/2,ab-0.2,0,"Rotation "+CStr(rot)
if flon=1 then
 dcCreateText (al+ar)/2+3,ab-0.2,0,"With Longitude correction"
else
 dcCreateText (al+ar)/2+3,ab-0.2,0,"Without Longitude correction"
end if
if feot=1 then
 dcCreateText (al+ar)/2-3,ab-0.2,0,"With EOT correction"
else
 dcCreateText (al+ar)/2-3,ab-0.2,0,"Without EOT correction"
end if

dcSetLineParms dcBlack, dcSolid, dcThin

dcSetTextParms dcBlack,"Tahoma","Bold",90,8,21,0,0
dcCreateText al-0.2,(at+ab)/2,0,"Author: Valentin Hristov, Sofia, Bulgaria      E-mail: valhrist@bas.bg"
dcSetTextParms dcBlack,"Tahoma","Bold",-90,8,21,0,0
dcCreateText ar+0.2,(at+ab)/2,0,"Web page: www.math.bas.bg/complan/valhrist/mystuff.htm"

'Central point of the hour lines
'dcSetCurrentLayer "ShadowCastingPoints"
'if abs(zNN0)>0.1 then
' dcCreateCircle -a*xNN0/zNN0,-b*yNN0/zNN0,.1
' dcCreateCircle -a*xNN0/zNN0,-b*yNN0/zNN0,.02
'end if

Double_Bifilar:

dcSetCurrentLayer "Gnomons"
dcSetCircleParms dcPurple,dcSolid,dcNormal
dcSetLineParms dcPurple,dcArrow,dcNormal
dccreateCircle g1px,g1py,.2
dcCreateLine g1px+g1vx,g1py+g1vy,g1px,g1py
dcSetCircleParms dcGreen,dcSolid,dcNormal
dcSetLineParms dcGreen,dcArrow,dcNormal
dccreateCircle g2px,g2py,.18
dcCreateLine g2px+g2vx,g2py+g2vy,g2px,g2py
dcSetCircleParms dcBlack,dcSolid,dcNormal
dcSetLineParms dcBlack,dcArrow,dcThin
dcSetCurrentLayer "default"
if per=0 then
 spring=1
 Drawing
end if
if per=1 then
 spring=0
 Drawing
end if
if per=2 then
 spring=1
 Drawing
 spring=0
 Drawing
end if
'if fdbl=1 then
 g1px=-g1px
 g1vx=-g1vx
 g2px=-g2px
 g2vx=-g2vx
if fdbl=1 then
 fdbl=0
 goto Double_Bifilar
end if

dcTurnLayerOff "DateLines"
dcTurnLayerOff "ShadowCastingPoints"
dcViewAll

End Sub 'Main

Sub Drawing
' 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 366 'max number of points in Spline is 248
' for count=1 to 366 step 2
 if spring=1 then nom=183 else nom=184
 for ind1=1 to nom
  if spring=1 then
   if ind1<12 then count=354+ind1 else count=ind1-11
  else
   count=171+ind1
  end if
  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
  
g1nx=g1vy*z4-g1vz*y4
g1ny=g1vz*x4-g1vx*z4
g1nz=g1vx*y4-g1vy*x4
l2=(g1nx*(g1px-g2px)+g1ny*(g1py-g2py)+g1nz*(g1pz-g2pz))/(g1nx*g2vx+g1ny*g2vy+g1nz*g2vz)
g2x=g2px+l2*g2vx
g2y=g2py+l2*g2vy
g2z=g2pz+l2*g2vz

'xx=g2x-g2z*x4/z4
'yy=g2y-g2z*y4/z4

g2nx=g2vy*z4-g2vz*y4
g2ny=g2vz*x4-g2vx*z4
g2nz=g2vx*y4-g2vy*x4
l1=(g1nx*(g2px-g1px)+g2ny*(g2py-g1py)+g2nz*(g2pz-g1pz))/(g2nx*g1vx+g2ny*g1vy+g2nz*g1vz)
g1x=g1px+l1*g1vx
g1y=g1py+l1*g1vy
g1z=g1pz+l1*g1vz

xx=g1x-g1z*x4/z4
yy=g1y-g1z*y4/z4

  if (z1<0 and z3<0 _
      and xx>al and xx<ar _
      and yy>ab and yy<at _
      and g2z>=0 and g1z>=0) then 'illuminated rectangle
   cs=cs+1
   spl(2*cs-1)=xx
   spl(2*cs)=yy

'''''''''''''''''''''''''''''''''
' All points on the gnomons which cast a shadow on the hour lines
' inside the rectangular area.
   dcSetCurrentLayer "ShadowCastingPoints"
   dcCreateCircle g1x,g1y,.05
   dcCreateCircle g2x,g2y,.05
   dcSetCurrentLayer "default"
'''''''''''''''''''''''''''''''''
   xe=spl(2*cs-1)
   ye=spl(2*cs)
   if hhb=True then
    xb=spl(2*cs-1)
    yb=spl(2*cs)
    if spring=1 then
     if ha=Int(ha) then
      dcSetTextParms dcBlue, "Tahoma","Bold",0,8,21,0,0
      dcCreateText xb,yb+.07,0,CStr(ham)
     end if
    else
     if ha=Int(ha) then
      dcSetTextParms dcRed, "Tahoma","Bold",0,8,21,0,0
      dcCreateText xb,yb-.07,0,CStr(ham1)
     end if
    end if 'spring
    hhb=False
    hhv=True
   end if
  else
   if cs>2 then
    if spring=1 then
     if ha=int(ha) then
      dcSetLineParms dcBlue, dcSOLID, dcTHICK
      dcSetSplineParms dcBlue, dcSOLID, dcTHICK
     else
      dcSetLineParms dcBlue, dcSOLID, dcTHIN
      dcSetSplineParms dcBlue, dcSOLID, dcTHIN
     end if
    else
     if ha=int(ha) then
      dcSetLineParms dcRed, dcSOLID, dcTHICK
      dcSetSplineParms dcRed, dcSOLID, dcTHICK
     else
      dcSetLineParms dcRed, dcSOLID, dcTHIN
      dcSetSplineParms dcRed, dcSOLID, dcTHIN
     end if
    end if
    dcCreateSpline spl(1),cs,False
    cs=0
   else
    if cs=2 then
     dcCreateLine spl(1),spl(2),spl(3),spl(4)
    'else
     'dcCreateCircle spl(1),spl(2),.02
    end if
    cs=0
   end if
  end if ' illuminated
 next ind1

 hhe=False
 if cs>2 then                       'if the hour line finishes at 365
  if spring=1 then
   if ha=int(ha) then
    dcSetLineParms dcBlue, dcSOLID, dcTHICK
    dcSetSplineParms dcBlue, dcSOLID, dcTHICK
   else
    dcSetLineParms dcBlue, dcSOLID, dcTHIN
    dcSetSplineParms dcBlue, dcSOLID, dcTHIN
   end if
  else
   if ha=int(ha) then
    dcSetLineParms dcRed, dcSOLID, dcTHICK
    dcSetSplineParms dcRed, dcSOLID, dcTHICK
   else
    dcSetLineParms dcRed, dcSOLID, dcTHIN
    dcSetSplineParms dcRed, dcSOLID, dcTHIN
   end if
  end if
  dcCreateSpline spl(1),cs,False
 else
  if cs=2 then
   dcCreateLine spl(1),spl(2),spl(3),spl(4)
  'else
   'dcCreateCircle spl(1),spl(2),.02
  end if
 end if
 if spring=1 then
  if ha=Int(ha) and hhe=False and hhv=True then
   dcSetTextParms dcRed, "Tahoma","Bold",0,8,21,0,0
   dcCreateText xe,ye-.07,0,CStr(ham1)
  end if
 else
  if ha=Int(ha) and hhe=False and hhv=True then
   dcSetTextParms dcBlue, "Tahoma","Bold",0,8,21,0,0
   dcCreateText xe,ye+.07,0,CStr(ham)
  end if
 end if
 cs=0

next ha

''''''''''''''''''''''''''''''''''''''''''''''''
dcSetCurrentLayer "DateLines"
''''''''''''''''''''''''''''''''''''''''''''''''
'Date lines
cs=0                         'counter for spline
if spring=1 then
 bb=1
 ee=19
else
 bb=19
 ee=37
end if
for i=bb to ee
 count=ind(i)
 hhv=False
 hhb=True
 hhe=False
 for ha=0 to 23.999 step .1
  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

g1nx=g1vy*z4-g1vz*y4
g1ny=g1vz*x4-g1vx*z4
g1nz=g1vx*y4-g1vy*x4
l2=(g1nx*(g1px-g2px)+g1ny*(g1py-g2py)+g1nz*(g1pz-g2pz))/(g1nx*g2vx+g1ny*g2vy+g1nz*g2vz)
g2x=g2px+l2*g2vx
g2y=g2py+l2*g2vy
g2z=g2pz+l2*g2vz

'xx=g2x-g2z*x4/z4
'yy=g2y-g2z*y4/z4

g2nx=g2vy*z4-g2vz*y4
g2ny=g2vz*x4-g2vx*z4
g2nz=g2vx*y4-g2vy*x4
l1=(g1nx*(g2px-g1px)+g2ny*(g2py-g1py)+g2nz*(g2pz-g1pz))/(g2nx*g1vx+g2ny*g1vy+g2nz*g1vz)
g1x=g1px+l1*g1vx
g1y=g1py+l1*g1vy
g1z=g1pz+l1*g1vz

xx=g1x-g1z*x4/z4
yy=g1y-g1z*y4/z4

  if (z1<0 and z3<0 _
      and xx>al and xx<ar _
      and yy>ab and yy<at _
      and g2z>=0 and g1z>=0) then 'illuminated rectangle
   cs=cs+1
   spl(2*cs-1)=xx
   spl(2*cs)=yy

'''''''''''''''''''''''''''''''''
' All points on the gnomons which cast a shadow on the hour lines
' inside the rectangular area.
   dcSetCurrentLayer "ShadowCastingPoints"
   dcCreateCircle g1x,g1y,.05
   dcCreateCircle g2x,g2y,.05
   dcSetCurrentLayer "DateLines"
'''''''''''''''''''''''''''''''''
    xe=spl(2*cs-1)
    ye=spl(2*cs)
    if hhb=True then
     xb=spl(2*cs-1)
     yb=spl(2*cs)
     if i+1=Int((i+1)/3)*3 then
      if spring=1 then
       dcSetTextParms dcBlue, "Tahoma","Bold",0,8,21,0,0
      else
       dcSetTextParms dcRed, "Tahoma","Bold",0,8,21,0,0
      end if
      dcCreateText xb,yb,0,"1 "+datetext((i+1)/3)
     end if
     hhb=False
     hhv=True
    end if
'   end if
  else
   if cs>2 then
      if spring=1 then
       dcSetTextParms dcBlue, "Tahoma","Bold",0,8,21,0,0
      else
       dcSetTextParms dcRed, "Tahoma","Bold",0,8,21,0,0
      end if
    if i+1=Int((i+1)/3)*3 then
     if spring=1 then
      dcSetLineParms dcBlue, dcSOLID, dcTHICK
      dcSetSplineParms dcBlue, dcSOLID, dcTHICK
     else
      dcSetLineParms dcRed, dcSOLID, dcTHICK
      dcSetSplineParms dcRed, dcSOLID, dcTHICK
     end if
     dcCreateText xe,ye,0,"1 "+datetext(int(i+1)/3)
    else
     if spring=1 then
      dcSetLineParms dcBlue, dcSOLID, dcTHIN
      dcSetSplineParms dcBlue, dcSOLID, dcTHIN
     else
      dcSetLineParms dcRed, dcSOLID, dcTHIN
      dcSetSplineParms dcRed, dcSOLID, dcTHIN
     end if
    end if
    dcCreateSpline spl(1),cs,False
    cs=0
   else
    if cs=2 then
     dcCreateLine spl(1),spl(2),spl(3),spl(4)
    'else
     'dcCreateCircle spl(1),spl(2),.02
    end if
    cs=0
   end if
  end if
 next ha

 hhe=False
 if cs>2 then
    if i+1=Int((i+1)/3)*3 then
     if spring=1 then
      dcSetLineParms dcBlue, dcSOLID, dcTHICK
      dcSetSplineParms dcBlue, dcSOLID, dcTHICK
     else
      dcSetLineParms dcRed, dcSOLID, dcTHICK
      dcSetSplineParms dcRed, dcSOLID, dcTHICK
     end if
    else
     if spring=1 then
      dcSetLineParms dcBlue, dcSOLID, dcTHIN
      dcSetSplineParms dcBlue, dcSOLID, dcTHIN
     else
      dcSetLineParms dcRed, dcSOLID, dcTHIN
      dcSetSplineParms dcRed, dcSOLID, dcTHIN
     end if
    end if
  dcCreateSpline spl(1),cs,False
 else
  if cs=2 then
   dcCreateLine spl(1),spl(2),spl(3),spl(4)
  'else
   'dcCreateCircle spl(1),spl(2),.02
  end if
 end if
 cs=0
' dcSetLineParms dcBlack,dcThin
next i

dcSetCurrentLayer "default"

End Sub 'Drawing

'''''''''''''''''''''''''''''''''''''''
Sub Latitude

dcCreateCircleEx -5,0,-5,2.5,-5,-2.5,2.5,2.5,0,0
dcCreateCircleEx 5,0,5,-2.5,5,2.5,2.5,2.5,0,0

for count=90 to 270 step 2
x=cos(count*d2r)
y=sin(count*d2r)
dcCreateLine -5+2.5*x,2.5*y,-5+2.55*x,2.55*y
next count

for count=90 to 270 step 10
x=cos(count*d2r)
y=sin(count*d2r)
dcCreateLine -5+2.5*x,2.5*y,-5+2.6*x,2.6*y
next count

for count=90 to 270 step 20
x=cos(count*d2r)
y=sin(count*d2r)
dcCreateLine -5+2.5*x,2.5*y,-5+2.65*x,2.65*y
next count

for count=-90 to 90 step 2
x=cos(count*d2r)
y=sin(count*d2r)
dcCreateLine 5+2.5*x,2.5*y,5+2.55*x,2.55*y
next count

for count=-90 to 90 step 10
x=cos(count*d2r)
y=sin(count*d2r)
dcCreateLine 5+2.5*x,2.5*y,5+2.6*x,2.6*y
next count

for count=-90 to 90 step 20
x=cos(count*d2r)
y=sin(count*d2r)
dcCreateLine 5+2.5*x,2.5*y,5+2.65*x,2.65*y
next count

if l>=0 then ll=l else ll=90+l

x=cos((90+2*ll)*d2r)
y=sin((90+2*ll)*d2r)
dcCreateLine -5+2.5*x,2.5*y,-5,2.5
dcCreateLine -5+2.5*x,2.5*y,-5,-2.5

x=cos((90-2*ll)*d2r)
y=sin((90-2*ll)*d2r)
dcCreateLine 5+2.5*x,2.5*y,5,2.5
dcCreateLine 5+2.5*x,2.5*y,5,-2.5

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