#   Portions Copyright (C) February 1998, Dave Fischer

#   dfischer@provide.net. SUN-PAL; V0.10; 2/1/98.

#   SUNCALC Port to CASL.  SunCalc was originally

#   written by Dave Rusin, rusin@math.niu.edu, on

#   4/10/94.  Suncalc and documentation	are found at: 

#   www.math.niu.edu/~rusin/uses-math/position.sun

#   it was written in UA Basic, a neat language which

#   is available at: ftp.coast.net/SimTel/msdos/ubasic/

#   This program, like SUNCALC, is freeware.  There is

#   NO WARRANTY and accuracy is not known.  Source code

#   may be distributed and used in other freeware

#   projects. Dave Rusin's UA Basic code was de-goto'd

#   and	made into functions, I keep his line number range

#   in the comments at the start of each function.

#   

#   Requires:

#   CASL runtime at: www.caslsoft.com

#   MathLib at: www.probe.net/~rhuebner/mathlib.html

#   Pilot OS 2.0 

#

#	Revisions:

#

#	V0.10 - 2/1/98  - Initial release.

#	

#	V1.00 - 3/20/98 - Fixed March equinox code.  Added

#                     Philadelphia to selector.



#   Variable definitions 



variables;



     Slant=23.45;       # angle of inclination of axis to earth orbit

     Latitude=42.22;    # data for Ann Arbor, MI, USA

     T_Lat = " ";

     Days=91;           # days since March equinox (assumed to be 3/21)

     DS = 0;			# DS=0 no daylight savings, DS=1 daylight savings

     T_Day = " ";

     Timex=12.0;        # hours since midnight at which you want to find the sun

     T_Time = " ";

     Nx=0;		        #(North,East,Up)

     Ex=0;		        #"up"=(0,0,1);"northwest"=(1,-1,0), etc.

     Ux=1.0;  

     Alpha=1.0;

     Theta=1.0;

     Psi=1.0;       

     Noon=1.0;

     Phix=1.0;

     pi=3.1415926535897;

     N=1.0;

     E=1.0;

     U=1.0;

     Ang1=1.0;

     Ang2=1.0;

     Ang3=1.0;

     Energy=1.0;

     Phi=1.0;

     X=1.0;

     A=1.0;

     B=1.0;

     C=1.0;

     F=1.0;

     Nroots=0;

     Phi1=1.0;

     Root1=1.0;

     Root2=1.0;

     Dawn=1.0;

     Dusk=1.0;

     nl=char(10);

     x_pos =30;

     y_pos= 30;

     flag1 = 0;

     flag2 = 0;

     degree = char(186);

     string bigstring;

     string str_arr[11] = "...";

     cons_indx = 0;

     sel_city[58] = "Ann Arbor-MI","Mobile-AL","Anchorage-AK","Phoenix-AZ","LittleRock-AR","LosAngeles-CA","SanDiego-CA","SanFrancisco-CA","Denver-CO","Hartford-CT","Wilmington-DE","Washington-DC","KeyWest-FL","Miami-FL","Pensacola-FL","Atlanta-GA","Honolulu-HI","Boise-ID","Chicago-IL","Indianapolis-IN","DesMoines-IA","Topeka-KS","Louisville-KY","NewOrleans-LA","PresqueIsle-ME","Baltimore-MD","Boston-MA","Detroit-MI","CopperHarb-MI","Minneapolis-MN","Jackson-MS","KansasCity-MO","Billings-MT","Lincoln-NE","LasVegas-NV","Manchester-NH","Newark-NJ","Albuquerque-NM","NewYork-NY","Charlotte-NC","GrandForks-ND","Columbus-OH","OklahomaCty-OK","Redmond-OR","Philadelphia-PA","Pittsburgh-PA","Providence-RI","Greenville-SC","SiouxFalls-SD","Memphis-TN","Dallas-TX","Houston-TX","SaltLakeCt-UT","Montpelier-VT","Richmond-VA","Seattle-WA","Wheeling-WV","Milwaukee-WI","Casper-WY";

     sel_lat[58] = 42.22,30.68,61.22,33.43,35.22,33.93,32.73,37.62,39.75,41.73,39.67,38.85,24.55,25.65,30.47,33.78,21.35,43.57,41.9,39.73,41.53,39.07,38.23,29.98,46.68,39.18,42.37,42.42,47.47,44.83,32.32,39.32,45.8,40.85,36.08,42.93,40.7,35.05,40.65,35.22,47.95,40,35.4,44.27,39.88,40.35,41.73,34.85,43.58,35.35,32.9,29.97,40.78,44.2,37.5,47.45,40.18,43.12,42.92;

     string str2;

     numeric ind1;

     x_org = 50;

     y_org = 130;

     x_up = 0;

     y_up = 0;

     sq_width = 252 - 63;

     sq_height = 315 - 63;

     circ_rad = 15;

     point_scale = 200;

     north = 0;

     east = 0;

     up = 0;

     ang2b = 0;

     old_north = 0;

     old_east = 0;

     old_up = 0;

     old_ang2 = 0;



end;

#    end variables           



#

#Main Frame "Working"

#

frame Working;

display "Sun-Pal";

end;



#

#Initial RUN button

#

button Execute, Working;

display "GO";

position  x_pos + 840, y_pos + 85;

pixel_size 100, 70;

end;



#

#Initial output "console"

#

text T1, Working;

position  x_pos, y_pos + 640;

char_size 17,3;

end;



#

#console UP (<) scroll button

#

button Cons_Up, Working;

display "<";

position  x_pos + 890, y_pos + 680;

pixel_size 60, 60;

end;



#

#console DOWN (<) scroll button

#

button Cons_Down, Working;

display ">";

position  x_pos + 890, y_pos + 760;

pixel_size 60, 60;

end;



#

# Selector1

#

# City Selector Object Definition

#

selector Selector1, Working;

	position x_pos + 465, y_pos + 470;

	pixel_size 480,150;

	list sel_city;

end;





#

#Latitude, defaults to 42.22 A2-MI

#

text T2, Working;

position  x_pos + 120, y_pos + 490;

char_size 4,1;

end;



label L2, Working;

position  x_pos,  y_pos + 490;

display "Lat: ";

end;



#

#Daylight savings button

#

button DaySave, Working;

display "DS";

position  x_pos + 720 , y_pos + 85;

pixel_size 100, 70;

end;



label L5, Working;

position  x_pos + 500,  y_pos + 170;

display "  ";

end;



#

#Get time button

#

button GetNow, Working;

display "Now";

position  x_pos + 785 , y_pos + 180;

pixel_size 160, 70;

end;



#

#About button

#

button About, Working;

display "About";

position  x_pos, y_pos + 870;

pixel_size 180, 70;

end;



#

#Time hours since midnight

#

text T3, Working;

position  x_pos + 650, y_pos + 340;

char_size 6,1;

end;



label L3, Working;

position  x_pos + 500,  y_pos + 340;

display "Time: ";

end;



#

#Days days since March Equinox (3/21, mostly)

#

text T4, Working;

position  x_pos + 650, y_pos + 260;

char_size 6,1;

end;



label L4, Working;

position  x_pos + 500,  y_pos + 260;

display "Days: ";

end;



#	 

#    SunCalc function 1. 630-720

#

function int_data_conv;

	 variables;

	     Length=1.0;

	 end;

	 Length=sqrt(Nx*Nx+Ex*Ex+Ux*Ux); 

     Nx=Nx/Length;

	 Ex=Ex/Length;

	 Ux=Ux/Length;

     Alpha=Slant*pi/180;

     Theta=Latitude*pi/180;

     Psi=Days*2*pi/365.25;

     Noon=atan(tan(Psi)*cos(Alpha));

     if cos(Theta)*cos(Psi)*cos(Noon)<0;

        Noon=Noon+pi;

	 end_if;

     Phix=Noon+(Timex-12)*2*pi/24;



end;

#   end function 1. 630-720



#

#    SunCalc function 2. 740-980

#

function app_solar_pos;

	variables;

	    N1=1.0;

		N2=1.0;

		N3=1.0;

		U1=1.0;

		U2=1.0;

		U3=1.0;

	 end;

     N1=-sin(Theta)*cos(Psi)*cos(Phi);

     N2=-cos(Alpha)*sin(Theta)*sin(Psi)*sin(Phi);

     N3=sin(Alpha)*cos(Theta)*sin(Psi);

     N=N1+N2+N3;

     E=-cos(Psi)*sin(Phi)+cos(Alpha)*sin(Psi)*cos(Phi);

	 U1=cos(Theta)*cos(Psi)*cos(Phi);

     U2=cos(Alpha)*cos(Theta)*sin(Psi)*sin(Phi);

     U3=sin(Alpha)*sin(Theta)*sin(Psi);

     U=U1+U2+U3;

     if E=0;

	    Ang1=pi/2;

	 end_if;

     if E<>0;

	    Ang3=atan(abs(N/E));

	 end_if;

	 if E>0;

        Ang1=atan(N/E);

	 end_if;

	 if E<0;

        Ang1=atan(N/E)+pi;

	 end_if;

     if U=1;

	    Ang2=pi/2;

     end_if;

	 if U=-1;

        Ang2=-pi/2;

	 end_if;

     if abs(U)<1;

	    Ang2=atan(U/sqrt(1-U*U));

	 end_if;

     if U<0;

	    Energy=0;

	 else;

	    Energy=(N*Nx+E*Ex+U*Ux)*100;

	 end_if;	

	 if Energy<0;

        Energy=0;

     end_if;



end;

#   end function 2. 740-980



#

#    SunCalc function 5. 1310-1380

#

function time_clean (numeric time_in) as string;



    variables;

	   numeric hour_part;

 	   numeric min_part;

    end;



    hour_part = int(time_in);

	time_in = time_in - hour_part;

	min_part = int(0.5+60*time_in);



	if min_part <10;

	   time_clean = string(hour_part,"") + ":0" + string(min_part,"");

	else;

       time_clean = string(hour_part,"") + ":" + string(min_part,"");

	end_if;



end;

#  end function 5. 1310-1380



#

#    SunCalc function 3. 1000-1150

#

function find_roots;



F=sqrt(A*A+B*B); 



if F=0;

   Nroots=0;

   flag2 = 1;

end_if;



if flag2 <> 1;

   A=A/F;

   B=B/F;

   C=-C/F;

end_if;



if B=0 and flag2 <> 1;

   X=A*pi/2;

else;

   X=atan(A/B);

   if B*cos(X)<0;

      X=X+pi;

   end_if;

end_if;



if abs(C)>1 and flag2 <> 1;

   Nroots=0;

   flag2 = 1;

end_if;



if abs(C)=1 and flag2 <> 1;

   Phi1=C*pi/2;

   Root1=Phi1-X;

   Nroots=1;

   flag2 = 1;

end_if;



if flag2 <> 1;

   Phi1=atan(C/sqrt(1-C*C));

   Root1=Phi1-X;

   Root2=pi-Phi1-X;

   Nroots=2;

end_if;



#reset flag2

flag2 = 0;



end;

#  end function 3. 1000-1150





#

# function draw_circ - draws a circle of size rad.

#

function draw_circ (numeric x_cent, numeric y_cent, numeric rad);

   set Working, pen, (x_cent-rad), (y_cent-rad);

   draw Working, circle, 2*rad, 2*rad;

end;



#

# function draw_angle - draws an angle

#

function draw_angle (numeric up, numeric method);



   variables;

      numeric vect;

   end;



   vect = sq_width - 80;

   set Working, pen, (x_org+sq_width+30), ((y_org+sq_height/2) - 10);



   x_up = vect*cos(abs(ang2b));

	

	if up>=0;

	   y_up = -vect*sin(abs(ang2b));

	else;

	   y_up = vect*sin(abs(ang2b));

	end_if;



   if method = 1;

      draw Working, line, x_up, y_up;

   else;

      clear Working, line, x_up, y_up;

   end_if;



end;



#

# function draw_sqr - draws a square

#

function draw_sqr(numeric x_up, numeric y_up, numeric sq_w, numeric sq_h);

   set Working, pen, x_up, y_up;

   draw Working, line, sq_w, 0;

   draw Working, line, 0, sq_h;

   draw Working, line, -sq_w, 0;

   draw Working, line, 0, -sq_h;

end;



#

# function draw_line method=1 draw; method=0 clear

# north changes to -north because +x drawing is

# lower on display



function draw_line(numeric method);



   north = - north * point_scale;

   east = east * point_scale;

   up = up * point_scale;



   if sqrt((north*north)+(east*east)) >= (sq_width/2);

       north = north/2;

       east = east/2;

   end_if;

   

   set Working, pen, x_org+(sq_width/2), y_org+(sq_height/2);

   

   if method = 1;

      draw Working, line, east, north;

      draw_angle (up, 1);

   else;

      clear Working, line, east, north;

      draw_angle (up, 0);

   end_if;



end;



#

# function draw_screen draws compass boxes

#

function draw_screen;

# draw two boxes

   draw_sqr(x_org, y_org, sq_width, sq_height);

   draw_sqr(x_org+sq_width+20, y_org, sq_width, sq_height);

   set Working, pen, (x_org+sq_width+20), (y_org+sq_height/2);

   draw Working, line, sq_width, 0;

# draw N for north

   set Working, pen, (x_org+sq_width/2-10), (y_org+20);

   draw Working, line, 0, -20;

   draw Working, line, 20, 20;

   draw Working, line, 0, -20;

# draw E for east

   set Working, pen, (x_org+sq_width), (y_org+sq_height/2-10);

   draw Working, line, -20, 0;

   draw Working, line, 0, 20;

   draw Working, line, 20, 0;

   set Working, pen, (x_org+sq_width-20), (y_org+sq_height/2);

   draw Working, line, 10, 0;



   draw_circ (x_org+(sq_width/2),y_org+(sq_height/2),circ_rad);

end;



#

# function Execute - the big, nasty working function

#

function Execute;



   variables;

       l = 0;

   end;



#  Get new T2, T3, T4 if user has changed

#  reset cons_indx, bigstring and str_arr[11]



   cons_indx = 0;



   while l <= 10;

         str_arr[l] = "...";

		 l = l + 1;

   end_while;



get T2, T_Lat;

Latitude = value(T_Lat);

get T3, T_Time;

Timex = value(T_Time);

get T4, T_Day;

Days = value(T_Day);



if DS = 1;

   Timex = Timex - 1.00;

end_if;



call int_data_conv;

str_arr[0] = "Latitude: " + string(Latitude, "0.##");

str_arr[1] = "Days since March Eqnx: " + string(Days, "");



#

#       show data for the moment 

#



if DS = 1;

   str_arr[2] = "At " + time_clean(Timex + 1.00) + " Sun at:";

else;

   str_arr[2] = "At " + time_clean(Timex) + " Sun at:";

end_if;



Phi=Phix;

call app_solar_pos;



str_arr[3] = "N=" + string(N,"0.###") + " E=" + string(E,"0.###") + " U=" + string(U,"0.###");



if U=1;

   str_arr[4] =	"overhead";

   flag1 = 1;

end_if;



if U=-1;

   str_arr[4] =	"underneath";

   flag1 = 1;

end_if;



if flag1 <> 1;

   str_arr[5] =	"Sighting is ";

   if E = 0 and N > 0;

      str_arr[5] = str_arr[5] +	"due N & ";

      str_arr[5] = str_arr[5] +	string(abs(Ang2)*(180/pi),"") + degree;

      if U>=0;

         str_arr[5] = str_arr[5] + " up.";

      end_if;

      if U<0;

         str_arr[5] = str_arr[5] + " down.";

      end_if;

   end_if;

   if E = 0 and N < 0;

      str_arr[5] = str_arr[5] +	"due S & ";

      str_arr[5] = str_arr[5] +	string(abs(Ang2)*(180/pi),"") + degree;

      if U>=0;

         str_arr[5] = str_arr[5] + " up.";

      end_if;

      if U<0;

         str_arr[5] = str_arr[5] + " down.";

      end_if;

   end_if;

   if E <> 0 and N > 0;

      str_arr[5] = str_arr[5] + "due ";

      str_arr[5] = str_arr[5] + string(abs(Ang3)*(180/pi),"") + degree + " N & ";

      str_arr[5] = str_arr[5] +	string(abs(Ang2)*(180/pi),"") + degree;

      if U>=0;

         str_arr[5] = str_arr[5] + " up.";

      end_if;

      if U<0;

         str_arr[5] = str_arr[5] + " down.";

      end_if;

   end_if;

   if E <> 0 and N < 0;

      str_arr[5] = str_arr[5] +	"due ";

      str_arr[5] = str_arr[5] + string(abs(Ang3)*(180/pi),"") + degree + " S & ";

      str_arr[5] = str_arr[5] + string(abs(Ang2)*(180/pi),"") + degree;

      if U>=0;

         str_arr[5] = str_arr[5] + " up.";

      end_if;

      if U<0;

         str_arr[5] = str_arr[5] + " down.";

      end_if;

   end_if;   

end_if;



#erase line routine	for direction boxes

north = old_north;

east = old_east;

up = old_up;

ang2b = old_ang2; 

call draw_line(0);



#draw line routine for direction boxes

north = N;

old_north = north;

east = E;

old_east = east;

up = U;

old_up = up;

ang2b = Ang2;

old_ang2 = ang2b;

call draw_line(1);



call draw_screen;



#reset flag1

flag1 = 0;

#

#     show data for the day 

#

A=cos(Theta)*cos(Psi);

B=cos(Alpha)*cos(Theta)*sin(Psi); 

C=sin(Alpha)*sin(Theta)*sin(Psi);

call find_roots;



if Nroots>0;

   flag1 = 1;

end_if;



if flag1 <> 1;

   str_arr[6] =	"Sun";

   if sin(Psi)*cos(Alpha-Theta)>0;

      str_arr[6] = str_arr[6] + " up all day.";

   else;

      str_arr[6] = str_arr[6] + " down all day.";

   end_if;

   flag2 = 1;

end_if;



# reset flag1

flag1 = 0;



if flag2 <> 1;

   Dawn=(Root1-Noon)*12/pi+12;

   Dawn=Dawn-24*int(Dawn/24);

   Dusk=(Root2-Noon)*12/pi+12;

   Dusk=Dusk-24*int(Dusk/24);



   if DS = 1;

      Dawn = Dawn + 1.00;

      Dusk = Dusk + 1.00;

   end_if;



   if Nroots = 1 or Nroots = 0;

      str_arr[6] = str_arr[6] + " at horizon at " + time_clean(Dawn) + ".";

      flag1 = 1;

   end_if;

 

   if flag1 <> 1;

      if Dawn > Dusk;

         X=Dawn;

         Dawn=Dusk;

         Dusk=X;

         Root1=Root2;

      end_if;

      str_arr[7] = "Day has " + string(Dusk-Dawn,"0.##") + " hrs of light."; 

      str_arr[8] = "Dawn at " + time_clean(Dawn) + ", dusk at " + time_clean(Dusk) + ".";

   end_if;



   Phi=Root1;

   call app_solar_pos;

   str_arr[9] = "Sunrise at bearing " + string(Ang1*180/pi,"0.##") + degree;  



   Phi=Noon;

   call app_solar_pos;

   str_arr[10] = "High sun at angle " + string(Ang2*180/pi,"0.##") + degree;



end_if;



bigstring = str_arr[cons_indx] + nl + str_arr[cons_indx + 1] + nl + str_arr[cons_indx + 2];

put T1, bigstring;



end;



#

# function Cons_Up - scrolls data in output

#                    text area UP.        

function Cons_Up;



   if cons_indx <=9 and cons_indx > 0;

      cons_indx = cons_indx - 1;

	  bigstring = str_arr[cons_indx] + nl + str_arr[cons_indx + 1] + nl + str_arr[cons_indx + 2];

      put T1, bigstring;

   end_if;



   if cons_indx = 10;

      cons_indx = cons_indx - 1;

	  bigstring = str_arr[cons_indx] + nl + str_arr[cons_indx + 1] + nl + "...";

      put T1, bigstring;

   end_if;

end;



#

# function Cons_Down - scrolls data in output

#                      text area DOWN.        

function Cons_Down;



   if cons_indx = 9;

      cons_indx = cons_indx + 1;

	  bigstring = str_arr[cons_indx] + nl + "..." + nl + "...";

      put T1, bigstring;

   end_if;



   if cons_indx = 8;

      cons_indx = cons_indx + 1;

	  bigstring = str_arr[cons_indx] + nl + str_arr[cons_indx + 1] + nl + "...";

      put T1, bigstring;

   end_if;



   if cons_indx <=7;

      cons_indx = cons_indx + 1;

	  bigstring = str_arr[cons_indx] + nl + str_arr[cons_indx + 1] + nl + str_arr[cons_indx + 2];

      put T1, bigstring;

   end_if;

end;



#

# function Selector1 - invoked when user taps a

#                    city from the selector1.        

function Selector1;



	get Selector1, ind1;

	put T2, string(sel_lat[ind1],"0.##");

end;



#

# function DaySave - invoked when user taps the

#                    daylight savings button.        

function DaySave;



    if DS = 0;

	   DS = 1;

	   put L5, "DaySave";

	else;

	   DS = 0;

	   put L5, "       ";

	end_if;

end;





#

#  returns 1 if the input year is a leap year

#  returns 0 if the input year is not a leap year

#  leapyear info: qa.pica.army.mil/~jferree/leapyear.html

#

function leap_year (numeric year_in) as numeric;



   if ((year_in/4) = int(year_in/4)) and ((year_in/100)	<> int(year_in/100)) or ((year_in/400) = int(year_in/400));

      leap_year = 1;

   else;

      leap_year = 0;

   end_if;

end;



#

#  returns fractional current time

#

function time_now as numeric;

#

    time_now = hour() + (minute()/60);

end;



#

#  returns the number of days since the March equinox

#

function day_mar_eqnx(numeric year_in, numeric month_in, numeric day_in) as numeric;



   variables;

       mths[12]=31,28,31,30,31,30,31,31,30,31,30,31;

   	   mar_eqnx = 21;   #assumed, some years it is 3/22

	   i=0;

       sum_of_days=0;

   end;



   if leap_year(year_in) = 1;

      mths[1] = 29;

   end_if;



   if month_in > 3;

      i=3;

      while i < (month_in - 1);

	    sum_of_days = mths[i] + sum_of_days;

		i=i+1;

	  end_while;

	  sum_of_days = sum_of_days + (mths[2]-mar_eqnx) + day_in;

   end_if;



   if month_in < 3;

      i=3;

      while i < 12;

	    sum_of_days = mths[i] + sum_of_days;

		i=i+1;

	  end_while;

	  

	  if month_in = 1;

	     sum_of_days = sum_of_days + (mths[2]-mar_eqnx) + day_in;

	  end_if;



	  if month_in = 2;

	     sum_of_days = sum_of_days + (mths[2]-mar_eqnx) + mths[0] +  day_in;

   	  end_if;

   end_if;



   if month_in = 3 and day_in < 21;

      i=3;

      while i < 12;

	    sum_of_days = mths[i] + sum_of_days;

		i=i+1;

	  end_while;

	  sum_of_days = sum_of_days + (mths[2]-mar_eqnx) + mths[0] + mths[1] + day_in;

   end_if;



   if month_in = 3 and day_in >= 21;

	  sum_of_days = day_in - mar_eqnx;

   end_if;



   day_mar_eqnx = sum_of_days;



end;



#

# function GetNow - invoked when user taps the

#                    Now button.        

function GetNow;

   put T4, string(day_mar_eqnx(year(), month(), day()), "###");

   put T3, string(time_now(),"0.##");

end;



#

# function About -   invoked when user taps the

#                    About button.        

function About;



   variables;

     l = 0;

   end;



   while l <= 10;

     str_arr[l] = "...";

     l = l + 1;

   end_while;

   

   cons_indx = 0;



   str_arr[0] =	"SUNCALC Port to CASL. SunCalc";

   str_arr[1] =	"was written by Dave Rusin";

   str_arr[2] =	"(rusin@math.niu.edu), it is";

   str_arr[3] =	"at: www.math.niu.edu/~rusin/";

   str_arr[4] =	"uses-math/position.sun/. This";

   str_arr[5] =	"CASL program is in the public";

   str_arr[6] =	"domain, no warranties + unknown";

   str_arr[7] =	"accuracy. More information at";

   str_arr[8] =	"www.provide.net/~dfischer/pilot";

   str_arr[9] =	"Questions: dfischer@provide.net";

   str_arr[10] = "SUN-PAL; V1.00; 3/20/98.";



   bigstring = str_arr[0] + nl + str_arr[1] + nl + str_arr[2];

   put T1, bigstring;

end;





#

# function startup - invoked when program starts

#                    default data set.        

function startup;

   put T2, string(Latitude, "0.##");

   put T3, string(Timex, "0.##");

   put T4, string(Days, "###");

   call draw_screen;

end;