'********************************************************
'* NOT FINISHED YET - VARIANT FOR TESTING               *
'* SDABIFIL.bas is a DeltaCad macro for producing a     *
'* Bifilar Sundial with two gnomons parallel to the     *
'* plate with Longitude Correction and EOT Correction   *
'* and ARBITRARY position of the plate, created by      *
'* Valentin Hristov (valhrist@bas.bg)                   *
'* 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 two gnomons are parallel to the plane and their projections are
'* indicated by the central vertical and horizontal lines.
'* Folding the upper, lower, left and right triangles will indicate the
'* heights of the gnomons (the distance to the plane).
'* With smaller heights more hour lines appear in the printed area.

'* You can use parallelipipeds with the suitable heights as gnomons and
'* then you will have to cut the printing along the central lines
'* and to adjust the four parts to the sides of the thick gnomons.

'* A second layer "ShadowCastingPoints" is added to show all points on
'* the gnomons which cast a shadow on the hour lines inside the
'* rectangular area.
'* Also the common point of the extensions of the hour lines (without
'* EOT correction) is shown in this layer.
'* Do not forget to switch OFF the layer before printing.

'* The design allows the sundial to have arbitrary position.

'* If the two heights of the gnomons are suitably choosen, then the hour
'* lines could be equiangular (15 degrees per hour). For Decl=0, Incl=0,
'* Rot=0, the ratio should be SIN(latitude).

'* In the particular case of equal heights 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 centre at the given height is enough.

'* If you want to have full analemmas as time lines, run the macro with
'* one of the periods and then run it again without closing DeltaCad
'* inside it with the other period. (15 min option is not suitable!)

'* Printing on A4 paper gives a really portable sundial.
'* You can use a copy machine to enlarge to A3 if you wish to have
'* a bigger size.

'* 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,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 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,200,226, "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,190,180,10, "Rectangular area: Bottom / Top / Left / Right"
 TextBox 15,200,37,10, .ab
 TextBox 60,200,37,10, .at
 TextBox 105,200,37,10, .al
 TextBox 150,200,37,10, .ar

 OKButton 42,212,37,12
 CANCELButton 116,212,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 = 0
prompt.rot = 0
prompt.per = 0
prompt.time_int = 1
prompt.eot_c = 0
prompt.lon_c = 0
prompt.a = 1
prompt.b = 0.676876
prompt.ab = -4
prompt.at = 4
prompt.al = -4
prompt.ar = 4


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 = prompt.l
lon = prompt.lon
cm = prompt.cm
d = prompt.d
i = prompt.i
rot = prompt.rot
a = prompt.a
b = prompt.b
at = prompt.at
ab = prompt.ab
al = prompt.al
ar = prompt.ar
per = prompt.per
time_int=prompt.time_int
eot_c=prompt.eot_c
lon_c=prompt.lon_c
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

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

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-(0-a),0,at
dcCreateLine 0,at-(0-a),0-a,at
dcCreateLine 0,at-(0-a),a,at

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

dcCreateLine ar-(0-b),0,ar,0
dcCreateLine ar-(0-b),0,ar,0-b
dcCreateLine ar-(0-b),0,ar,b

dcSetCurrentLayer "default"

dcCreateCircle 0,ab-a,.05
dcCreateCircle 0,at-(0-a),.05
dcCreateCircle al-b,0,.05
dcCreateCircle ar-(0-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,10,21,0,0
dcCreateText (al-(0-ar))/2,at+0.2,0,p
dcCreateText (al-(0-ar))/2-2.5,at+0.2,0,"Latitude "+ll
dcCreateText (al-(0-ar))/2+2.5,at+0.2,0,"Longitude "+lol
dcCreateText (al-(0-ar))/2+1.5,ab-0.2,0,"Declination "+CStr(d)
dcCreateText (al-(0-ar))/2-1.5,ab-0.2,0,"Inclination "+CStr(i)
dcCreateText (al-(0-ar))/2,ab-0.2,0,"Rotation "+CStr(rot)
if flon=1 then
 dcCreateText (al-(0-ar))/2+3,ab-0.2,0,"With Longitude correction"
else
 dcCreateText (al-(0-ar))/2+3,ab-0.2,0,"Without Longitude correction"
end if
if feot=1 then
 dcCreateText (al-(0-ar))/2-3,ab-0.2,0,"With EOT correction"
else
 dcCreateText (al-(0-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-(0-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-(0-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 (0-a)*xNN0/zNN0,(0-b)*yNN0/zNN0,.1
 dcCreateCircle (0-a)*xNN0/zNN0,(0-b)*yNN0/zNN0,.02
end if
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

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
  
  if (z1<0 and z3<0 and (0-a)*x4/z4>0-(0-al) and (0-a)*x4/z4<ar and (0-b)*y4/z4>0-(0-ab) and (0-b)*y4/z4<at) then 'illuminated
   cs=cs+1
   xx=(0-a)*x4/z4
   yy=(0-b)*y4/z4
   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"
   dcCreateLine -.02,(a-b)*y4/z4,.02,(a-b)*y4/z4
   dcCreateLine (b-a)*x4/z4,-.02,(b-a)*x4/z4,.02
   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

''''''''''''''''''''''''''''''''''''''''''''''''
'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
  
  if (z1<0 and z3<0 and (0-a)*x4/z4>0-(0-al) and (0-a)*x4/z4<ar and (0-b)*y4/z4>0-(0-ab) and (0-b)*y4/z4<at) then 'illuminated
   cs=cs+1
   xx=(0-a)*x4/z4
   yy=(0-b)*y4/z4
   spl(2*cs-1)=xx
   spl(2*cs)=yy
    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
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