'********************************************************
'* NOT FINISHED YET - VARIANT FOR TESTING               *
'* CYLARB1.bas is a DeltaCad macro for producing a      *
'* Cylinder Sundial with arbitrary orientation and a    *
'* hole as a nodus and hour lines on the inner surface  *
'* with Longitude Correction and EOT Correction         *
'* Created by Valentin Hristov (valhrist@bas.bg) in the *
'* end of January 2008.                                 *
'* (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.                         *
'********************************************************

'* This Cylinder Sundial is of nodus type. The light goes through a hole
'* and the time can be found on the inner surface of the cylinder.
'* The rectangle on the drawing represents the whole surface. You have
'* to curve it until the upper and the lower boundaries coincide. On
'* this common line is the place of the hole (nodus).

'* The cylinder can have arbitrary orientation. To give the parameters,
'* imagine that it lies on a plane in horisontal left-right position,
'* The nodus hole is opposite to the plane. Then giving suitable
'* declination, inclination, and rotation of the plane, you can have
'* arbitrary orientation of the cylinder.

'* Interesting particular cases (northern hemisphere):
'*      decl    incl    rot
'*      0       lat     90      polar cylinder, hole to celestial South
'*      90      0       -lat    the same
'*      -90     0       lat     the same
'*      0       lat     0       E-W horizontal cylinder, hole to celestial South
'*      0       0       0       E-W horiz. cylinder, hole to Zenith
'*      180     0       0       the same with left<->right change
'*      0       90      0       E-W horiz. cylinder, hole to South
'*      0       90      90      vertical cylinder, hole to South
'*      90      90      90      vertical cylinder, hole to West
'*      90      90      0       N-S horiz. cylinder, hole to West
'*      90      0       0       N-S horiz. cylinder, hole to Zenith
'*      -90     0       0       the same with left<->right change
'*      -90     90      0       N-S horiz. cylinder, hole to East
'*      -90     90      90      vertical cylinder, hole to East
'* etc.

'* Different layers are drawn when you run the macro. To switch them ON
'* or OFF go to the "View-Layers" menu.

'* Before printing switch OFF the not needed layers.

'* 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,x5,y5,z5,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
call Sundial
'**********
Sub Sundial
init_constants
Input_constants_of_sundial
if button=0 then goto cncl
Angles
Vectors
Preparation
Main
'Latitude
cncl:
End Sub
'End of program


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

Begin Dialog CONSTANTS_INPUT 10,25,205,206, "Input data for the sundial"
 Text 15,0,180,10, "CYLINDER SUNDIAL WITH ARBITRARY ORIENTATION"
 Text 15,10,180,10, "with corrections for the longitude 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, "Left end"
 Text 115,170,180,10, "Right end"
' TextBox 15,180,37,10, .ab
' TextBox 60,180,37,10, .at
 TextBox 50,170,37,10, .al
 TextBox 150,170,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)"

 Text 55,180,150,10, "Cylinder radius"
 TextBox 110,180,37,10, .r

 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 = -8
'prompt.at = 8
prompt.al = -4
prompt.ar = 4
'prompt.g1px = 0
'prompt.g1py = 0
'prompt.g1pz = 2
'prompt.g1vx = 0
'prompt.g1vy = 1
'prompt.g1vz = 0
'prompt.g2px = 0
'prompt.g2py = 0
'prompt.g2pz = 2
'prompt.g2vx = 1
'prompt.g2vy = 0
'prompt.g2vz = 0
'prompt.dbl = 1
prompt.r = 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
r=CDbl(prompt.r)
at=R*pi
ab=-R*pi
g1px=0
g1py=0
g1pz=2*r
g1vx=0
g1vy=1
g1vz=0
g2px=0
g2py=0
g2pz=2*r
g2vx=1
g2vy=0
g2vz=0

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

dcCreateCircle 0,at,.1
dcCreateCircle 0,ab,.1

'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"
'dcTurnLayerOff "Gnomons"
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
xx=-2*r*x4*z4/(y4*y4+z4*z4)
yy=2*r*Atn(-y4/z4)

'  cylinder

  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
xx=-2*r*x4*z4/(y4*y4+z4*z4)
yy=2*r*Atn(-y4/z4)

'  cylinder


  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