unit U_SunPos2;
{Copyright 2001, Gary Darby, Intellitech Systems Inc., www.DelphiForFun.org

 This program may be used or modified for any non-commercial purpose
 so long as this original notice remains in place.
 All other rights are reserved
 }

{SunPos displays solar position and sunrise/sunset infor for any given date
 and time at any given location on earth.
 Also displays a plot of the shadow analemma for a location and time of day
 }

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls, NumEdit, ComCtrls, ExtCtrls;

type

  T3DVector=record
    x,y,z:extended;
end;

  TDrawMode=(sunlines, analemma, none );

  TForm1 = class(TForm)
    ShowBtn: TButton;
    Memo1: TMemo;
    AnalemmaBtn: TButton;
    Panel1: TPanel;
    Label2: TLabel;
    LongEdt: TEdit;
    EWRGrp: TRadioGroup;
    NSRGrp: TRadioGroup;
    LatEdt: TEdit;
    Label1: TLabel;
    Label3: TLabel;
    DatePicker: TDateTimePicker;
    Label8: TLabel;
    TZBox: TComboBox;
    Label9: TLabel;
    TimePicker: TDateTimePicker;
    DLSRGrp: TRadioGroup;
    PBox: TPaintBox;
    AnTypegrp: TRadioGroup;
    procedure ShowBtnClick(Sender: TObject);
    procedure FormActivate(Sender: TObject);
    procedure AnalemmaBtnClick(Sender: TObject);
    procedure DatePickerUserInput(Sender: TObject;
      const UserString: String; var DateAndTime: TDateTime;
      var AllowChange: Boolean);
    procedure PBoxPaint(Sender: TObject);
    procedure LongEdtExit(Sender: TObject);
  private
{ Private declarations }
public
{ Public declarations }

Adate:TDateTime;
    JDate:integer;
    tzHours:Integer; {time zone add to Greewich time. west = neg}
dlshours:integer;  {Daylight savings - hours added to Local time}
LocalTime:TDateTime;
    drawmode:TDrawMode;
    long,lat:extended;

{Variables used or set by SunPos function}

{Sunpos outputs}
sunrise:single; { Sunrise time                 }
sunset:single; { Sunset time                  }
sunriseangle, sunsetangle:extended;
    az,alt:extended;
    midday,middayaz,middayalt:extended;
    errmsg:string;
    Hourangle:extended;
    declination:extended;     { Angle of declination         }
EqOfTimeMinutes:extended;

{Drawing variables}
cx,cy,clen:integer;
    anpoints:array of tPoint;  {analemma points for drawing}
anscale:single; {analemma drawing scaling factor}
anplotx, anploty:integer; {coordinates for plot of "date" on analemma}

function sunpos( longitude, latitude, hrtime:extended;
                 tz,day:integer; var pos:T3DVector):boolean;
    procedure drawsunlines;
    function GetBaseData:boolean;
end;


var
  Form1: TForm1;

implementation

{$R *.DFM}
Uses math;

{************************** Local Routines ******************}

{degree trig functions}
function cosd(t:extended):extended;
  begin  result:= cos(0.01745329251994 * t); end;

  function sind(t:extended):extended;
  begin  result:= sin(0.01745329251994 * t); end;

  function tand(t:extended):extended;
  var  x:extended;
  begin
    x:=degtorad(t);
    result:= tan(x);
end;



{Angle/time  to/from string}

{**************** AngleToStr ***************}
function AngleToStr(angle:extended):string;
       var
         D:integer;
         M,S:extended;
       begin
         d:=Trunc(angle);
         m:=abs(frac(angle)*60);
         s:=frac(M)*60;
         m:=int(M);
         if s>=60-1e-10 then
         begin
           s:=s-60;
           m:=m+1;
end;
         if M=60 then
         begin
           m:=0;
           if d>=0 then inc(d) else dec(d);
end;
         result:=format('%3d %2d %5.2f',[d,trunc(m),s]);
end;

{****************** StrToAngle ******************}
function StrToAngle(s:string; var angle:extended):boolean;
       var
         n, sign:integer;
         w,ds,ms,ss:string;
       begin
         result:=true;
         ds:='0';
         ms:='0';
         ss:='0.0';
         w:=trim(s);
         if w[1]='-' then
         begin
           sign:=-1;
           delete(w,1,1);
end
else sign:=1;
         n:=pos(' ',w);
         if n>0 then
         begin
           ds:=copy(w,1,n-1);
           delete(w,1,n);
           w:=trim(w);
           n:=pos(' ',w);
           if n>0 then
           begin
             ms:=copy(w,1,n-1);
             delete(w,1,n);
             ss:=trim(w);
end
else ms:=w;
end
else ds:=w;
         try
           angle:=sign*(strtoint(ds)+strtoint(ms)/60+strtofloat(ss)/3600);
         except
           result:=false;
           angle:=0;
end;
end;

{Day of year function}
function JulianDay(d:TDateTime):integer;
      var
        y,m,day:word;
      begin
        decodedate(d,y,m,day);
        result:=trunc(d-encodedate(y-1,12,31));
end;


{****************************   SUNPOS *****************************}
function TForm1.sunpos( longitude, latitude, hrtime:extended;
                 tz,day:integer;
                 var pos:T3DVector):boolean;

  (*  Adapted from C code written by :
  /*                                                                      */
  /*                Joe Cychosz                                           */
  /*                Purdue University CADLAB                              */
  /*                Potter Engineering Center                             */
  /*                W. Lafayette, IN  47907                               */
  /*                                                                      */
  /* -------------------------------------------------------------------- */
  *)

     function eqoftime(jday:integer):extended;
{Equations of time - hours that solar differs from mean tine an any day}
var b:extended;
     begin
        b:= 360.0 * (day-81) / 364;
        result:=(9.87*sind(2.0*b) - 7.53*cosd(b) - 1.5*sind(b)) / 60.;  {eq of time}
end;

     (*
     function eqOfTime2(jday:integer):extended;
   var
     D,T,R0,R1:Extended;
     SG,TL,t0,UT: Extended;
     Flag:integer;
   begin
     D:=JDay;
     T:=D/365.25 -1;
     R0:=T*(5.13366E-02+T*(2.586222E-5-T*1.722E-9));
     R1:=6.697374558+2400.0*(T-(year-2000)/100.0);
     T0:=Adjust24(R0+R1, Flag);
     UT:=adjust24(value*24-FDLSHours-FTZHours, Flag);
     FUniversalTime:=UT/24.0+flag;
     T:=T0;
     If Flag>0 then T:=T+6.57098244E-2
     else if flag<0 then T:=T-6.57098244E-2;
     SG:=adjust24(UT*1.002737908+T, Flag);  {Greenwich Sidereal}
result:=Adjust24(sg+Fgeolonlat.x/15.0, Flag); {Local Sidereal}
end;
  *)


 var
    lontz:extended;   { time zone longitude}
b:extended;       { Time angle                   }
e:extended;       { Time                         }
{Ha:extended; }     { Hour angle                   }
we,gs:extended;
    phi:extended;     {altitude }
cws,ws:extended;      { Sunrise, sunset angle        }
x:extended;
{x1,x2,x3,x4,x5,x6:extended;}
temppos:t3dvector;
begin
   result:=false;
   errmsg:='';
   lontz:=tz*15; {time zone in degrees}
e:=eqoftime(day);
   midday:=12-long/15+tzhours-e{+dlsHours-e};
   while midday<0 do midday:=midday+24.0;
{sunpos(longitude, latitude, midday, tz,day,temppos);}
declination := 23.45 * sind (360.0 * (284+day) / 365.25);
   Hourangle:=180.0 + lontz - longitude - 15.0*(midday+e);
   phi:=arccos(sind(latitude)*sind(declination)
            + cosd(latitude)*cosd(declination)*cosd(Hourangle));{alt}
gs:=arcsin(cosd(declination) * sind(HourAngle) / sin(phi)); {az}
if  ( lat-declination)<0 then gs :=  pi - gs ;
   middayaz:=180-radtodeg(gs);
   middayalt:=90-radtodeg(phi);
   midday:=midday+dlshours;
{Compute solar declination (23.45 deg on June 21) }
declination := 23.45 * sind (360.0 * (284+day) / 365.25);
{Compute the hour angle in degrees. }
Hourangle:=180.0 + lontz - longitude - 15.0*(hrtime-dlshours+e);
   phi:=arccos(sind(latitude)*sind(declination)
              + cosd(latitude)*cosd(declination)*cosd(Hourangle));{alt}

x:=  -tand(Declination)*tand(Latitude) ; {replacement}
cws:=x+cosd(90.833)/(cosd(latitude)*cosd(declination));
   gs:=arcsin(cosd(declination) * sind(HourAngle) / sin(phi)); {az}
if  ( lat-declination)<0 then gs :=  pi - gs ;
   az:=180-radtodeg(gs);
   alt:=90-radtodeg(phi);
{ Compute the position of the sun. }
pos.x :=  sin (gs);    {East}
pos.y :=   -cos (gs);  {North}
pos.z :=   cos (phi);                    {Up}
EqOfTimeMinutes:=e*60;

   if cws<-1 then errmsg:='Sun never sets'
else if cws>1 then errmsg:='Sun never rises'
else
   begin
     result:=true;
{Compute sunrise, sunset times and angles}
ws:=radtodeg(arccos(cws));

     x:=(180.0 - ws + lontz - longitude);
     sunrise:=  x/ 15.0 - e+ DLSHours;
       x:=(180.0 + ws + lontz - longitude);
       sunset:=  x/ 15.0 - e + DLSHours;
{ result:=(hrtime >= sunrise) and (hrtime <= sunset);}
phi:=arccos(sind(declination)/cosd(latitude));
       sunriseangle:=radtodeg(phi);
       sunsetangle:=360-sunriseangle;
end;
end;

{**************** FormActivate ****************}
procedure TForm1.FormActivate(Sender: TObject);
begin
  TZBox.itemindex:=7;
  datepicker.date:=date;
  drawmode:=none;
  doublebuffered:=true;
  cx:=pbox.width div 2;
  cy:=pbox.height div 2;
  clen:=3*pbox.width div 8; {lines 3/4 out from center}
end;

{******************** DatePickerUserInput  *********}
procedure TForm1.DatePickerUserInput(Sender: TObject;
  const UserString: String; var DateAndTime: TDateTime;
  var AllowChange: Boolean);
{to overcome some default date entry edit problems}
begin
  try
    DateAndTime:=strtodate(userstring);
  except
end;
end;

{********************** GetBaseData ***************}
function TForm1.GetBaseData:boolean;
begin
  begin
    result:=true;
    Adate:=datepicker.date;
    JDate:=Julianday(Adate);
    tzHours:=tzBox.itemindex-12;
    if strtoangle(longedt.text,long) then
    begin
      Long:=abs(Long);
      if EWRGrp.itemindex>0 then Long:=-Long;
end
else
    begin
       showmessage('Invalid longitude, must be 1 to 3 numbers separated by spaces');
       result:=false;
end;
    if strtoangle(latedt.text,Lat) then
    begin
      Lat:=abs(Lat);
      if NSRGrp.itemindex>0 then Lat:=-Lat;
end
else
    begin
      showmessage('Invalid latitude, must be 1 to 3 numbers separated by spaces');
      result:=false;
end;
    dlshours:=DLSRGrp.itemindex;
    LocalTime:=(timepicker.time);
end;
end;

{************************** ShowBtnClick ********************}
procedure TForm1.ShowBtnClick(Sender: TObject);
{ Angles are in degrees, times are in hours. }
var
{Direction cosines of a vector pointing toward the sun from the origin.}
pos:T3dVector;
begin
   If GetBaseData then
   begin
      if drawmode=analemma then memo1.lines.clear;
      sunpos(long,lat,Timepicker.time*24,
            tzHours, Julianday(datepicker.date), pos);
      if errmsg='' then drawsunlines;
      with memo1,lines do
      begin
        add('**********************************');
        add('SUNPOS: Longitude:'+angletostr(long)
                    +', Latitude:' + angletostr(Lat)
                    + '  '+datetostr(datepicker.date)
                    );
        if errmsg='' then
        begin
          add('Sunrise:     '
+formatdatetime('hh:nn AM/PM',sunrise/24)
                 +format(', Azimuth %6.1f degrees from North',[sunriseangle]));
          add('Sunset:      '
+formatdatetime('hh:nn AM/PM',sunset/24)
                 +format(', Azimuth %6.1f degrees from North',[sunsetangle]));
end
else  add('For '+datetostr(datepicker.date)+' the sun '+errmsg);
        add('Solar noon: ' + timetostr(midday/24)+', Azimuth: '
+angletostr(middayaz)+', Altitude:'+angletostr(middayalt));

        add(format('Eq of time:  %4.2f minutes',[EqOfTimeMinutes])
                   + ', Declination: ' + angletostr(declination)+' deg');
       
        add('');
        add(timetostr(timepicker.time)
            + format(' local time,  Azimuth: %5.1f, Altitude: %5.1f,  ',[az,alt]));
        add(format('Position Vector (E,N,Up): %6.3f,%6.31f,%6.3f', [pos.x, pos.y, pos.z]));
end;
end;
end;

{********************* AnalemmaBtnClick **************************}
procedure TForm1.AnalemmaBtnClick(Sender: TObject);
{calculate the analemma as a shadow cast by a vertical pole
 plot top = south}
var
  i,n,nx,ny:integer;
  dayinc:integer;
  pos:T3DVector;
{L:extended; }
Hour:extended;
  xmax, ymax, xmin, ymin:integer;
  range:integer;
  nbrpoints:integer;
  anscaley,anscalex:extended;
  adjustx, adjusty:integer;
  pointto:extended;

        procedure makepoint(az,alt:extended; var x,y:integer);
          var
            L:extended;
          begin
          case antypegrp.itemindex of
            0: {shadow}
begin
              L:=1000/tand(alt); {defines a circle (x^2+y^2=L^2) upon which lies the point}
x:=-trunc(L*sind(az-180));
              y:=trunc(L*cosd(az-180));
end;
            1: {camera south}
begin
              x:=trunc(1000*sind(180-az));
              y:=trunc(1000*cosd(alt));
end;
            2: {camera at sun on input date time}
begin
              x:=trunc(1000*sind(pointto-az));
              y:=trunc(1000*cosd(alt));
end;
end; {case}
end;

begin
  if getbasedata then
  Begin
    with memo1.lines do
    begin
      clear;
      add('An analemma describes the shape of a figure formed by the sun''s');
      add('position in the sky if observed at the same time each day for a year.');
      add('');
      add('A common way to do this is to plot the tip of the shadow cast on the');
      add('ground by vertical stick.  This plot simulates that shadow as well as');
      add('a camera view with camera fixed pointing south or rotated to sun postition.');
      add('Camera fixed south produces "fisheye" lens effect for early or late observation times');
      add('');
      add('The analemma looks best when plotted for hours around noon.  It is');
      add('automatically scaled and panned to fit in the plot area but when the sun gets near');
      add('the horizon, shadow lengths become very large.  This can result in');
      add('truncated plot results');
      add('');
      add('The Red dot on analemma reflects the sun''s position on the date and time');
      add('entered in the input area');
end;
    drawmode:=analemma;
    dayinc:=5; {sample interval in days}
hour:=frac(localtime)*24; {which hour of the day to plot}
nbrpoints:= 365 div dayinc+1;
    setlength(anpoints, nbrpoints);
    xmax:=0; ymax:=0;
    xmin:=maxint; ymin:=maxint;
    sunpos(long, Lat,hour, trunc(tzhours),  JDate, pos);
    pointto:=az;  {direction of camera}
{plotindex:= trunc(0.0+round(JDate / dayinc));}
n:=0;
    i:=0;
    while i<nbrpoints do
    begin
      inc(i);

         sunpos(long, Lat,hour, trunc(tzhours),  dayinc*i, pos);
{only generate points if altiude > 0
         - before shadow tries to go to infinite length}
if alt>0 then
      begin
        inc(n);

        with anpoints[n-1] do
        begin
          makepoint(az,alt,x,y);
          xmax:=max(x,xmax);
          ymax:=max(y,ymax);
          xmin:=min(x,xmin);
          ymin:=min(y,ymin);
end;
end;
end;
    setlength(anpoints,n);
{set scale}
range:=ymax-ymin;
    anscaley:=0.9*pbox.height/range;
    range:=xmax-xmin;
    anscalex:=0.9*pbox.width/range;
    anscale:=min(anscalex,anscaley);
    if anscale<0.01 then anscale:=0.01;

    adjustx:=(xmax+xmin) div 2;
    adjusty:=(ymax+ymin) div 2;
{and scale the values for plotting}
for i:= 0 to high(anpoints) do
    with anpoints[i] do
    begin
      x:=cx+trunc(anscale*(x-adjustx));
      y:=cy+trunc(anscale*(y-adjusty));
end;
{calc plot coordinates for date entered by user}
sunpos(long, Lat,hour, trunc(tzhours),  JDate, pos);
    makepoint(az,alt,nx,ny);
    anplotx:=cx+trunc(anscale*(nx-adjustx));
    anploty:=cy+trunc(anscale*(ny-adjusty));
end;
  pbox.invalidate;
end;

{******************** DrawSumLines *************}
procedure TForm1.drawsunlines;
 begin
    drawmode:=sunlines;
    pbox.invalidate;
end;

{****************** PBoxPaint ****************}
procedure TForm1.PBoxPaint(Sender: TObject);
var
  i:integer;
  Rad,C:integer;
  CosA1,SinA1,CosA2,SinA2:extended;
  sign:Integer;
  begin
  if drawmode=none then exit;
  with pbox.canvas do
  begin
    If drawmode=sunlines then
    begin
      Rad:=70;
      IF round(middayaz)=180 then sign:=+1 else sign:=-1;
      brush.color:=clOlive;
      pen.color:=clblack;
      rectangle(rect(0,0,width,height));
      ellipse(cx-Rad, cy-Rad, cx+Rad, cy+Rad);
      font.size:=12;
      font.style:=[fsbold];
      textout(cx-10, cy-Rad-25,'N');
      textout(cx-10, cy+Rad,'S');
      textout(cx-Rad-25, cy-10,'W');
      textout(cx+Rad+10, cy-10,'E');
      ellipse(cx-2,cy-2,cx+2,cy+2);
      moveto(cx,cy);
      CosA1:=cos(degtorad(sunriseangle-90));
      SinA1:=sin(degtorad(sunriseangle-90));
      CosA2:=cos(degtorad(sunsetangle-90));
      SinA2:=sin(degtorad(sunsetangle-90));
      lineto(trunc(cx+clen*CosA1),
              trunc(cy+clen*SinA1));
      moveto(cx,cy);
      lineto(trunc(cx+clen*cosA2),
              trunc(cy+clen*sinA2));
      (*  attempt to plot ellipse sun area with tip proportional to altitude,
          not working for southern hemisphere
      c:=Rad*trunc(middayalt) div 90;
      arc(cx-trunc(Rad*cosA1), cy-sign*c+2*sign*trunc(rad*sinA1),
          cx+trunc(Rad*cosA1), cy+sign*c,
          cx-trunc(rad*cosA1), cy+trunc(rad*sinA1),
          cx+trunc(rad*cosA1), cy+trunc(Rad*sinA2));
       *)
       brush.color:=clyellow;
       floodfill(cx,cy+5*sign,clblack,fsborder);
end
else if drawmode=analemma then
    begin
      brush.color:=clblue;
      pen.color:=clblack;
      rectangle(rect(0,0,width,height));
      pen.color:=clyellow;
      with anpoints[high(anpoints)] do moveto(x,y);
      for i:=low (anpoints) to high(anpoints) do
      with anpoints[i] do lineto(x,y);
{draw current sun position}
pen.color:=clred;
      ellipse(anplotx-2,anploty-2,anplotx+2,anploty+2);
end;
end;
end;

procedure TForm1.LongEdtExit(Sender: TObject);
{Autoset the time zone when longitude changes}
var W:extended;

begin
  if strtoangle(longedt.text,W) then
  begin
    if ewrgrp.itemindex=1 then W:=-w;
    If (W<>long) then tzbox.itemindex:=12+trunc(round(w) / 15); {guess 15 degrees per time zone}
end;
end;


end.