{~ Planar/Cylindric Sundials and Insolation of Collectors and Windows ~}

program SunDial;   {by Jan Hollan, N.Copernicus Obs. and Planetarium, Brno}
(*  Copyright (C) 1999 Jan Hollan

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *)
uses solar_ut,angles_o,graph_m8,Dos,crt,g_scree,Str_Num,
     Un_Str_M {instead of unit_str by Blaise Computing, Inc.},
     params_g;

const mode:byte=9;
h1=cl+
'SunDial h  A#  h#  F#  DL#  X# XS# Y# YS# TA#  R#  {S/B} W[output file name]'+cl+
'   L# Ddd[.mm[.yy] Thh[.mm[.ss]] UT# AL# P# {Left/Right} M#'+cl+
'   E[<name>[:<name>]] AB# N# Q'+cl+
'  { PT[LaTeX picture file name] / PP[Encapsulated PostScript file name] }';
h2=
'plots a sundial on the screen. Parameters (on the first line) are:'+cl2;

h2a=
'H   writes this information into SunDial.hlp'+cl+
'A   Azimuth of the direction perpendicular to the sundial plane (default 180=S)'+cl+
'h#  angular height of that direction /1 degree (default 0, i.e. vertical wall)'+cl;
h3=
'F#  latitude (in degrees, too; default '
 +MyLocation_latitude_string+' = latitude of '
  +MyLocation_Name+')'+cl+
'DL# How many minutes is the local time ahead of that which you want to plot'+cl+
'X#  horizontal shift of the plot in percent of the screen width'+cl;
h35=
'XS# side along X coordinate (default 100 per cent of the screen)'+cl+
'Y#  vertical   shift of the plot in percent of the screen height'+cl+
'YS# side along Y coordinate (def. 94 % for VGA to match A4 landscape)'+cl;
h4=
'TA# is the "hour" along which the analema (the "8") is to be plotted'+cl+
'R#  is the "Radius of projection" = magnification (default 1)'+cl+
'C#  implies a cylinder instead of a plane; # is its radius relative to R'+cl;
h4a=
'S   asks for a Short description instead of a long one.'+cl;
h5=
'B   asks for a short description on the Bottom of the screen'+cl+
'BR- plots no border (BC would plot a circular border)'+cl+
'W<name> asks for making an output file <name> (default SunDial.txt)'+cl;

h6=
'The parameters on the second line concern the case, that the Azimuth'+cl+
'of the plane is not known, but its angle towards the Sun in some'+cl+
'instant is. The time is to be given in that case at least; the other parameters'+cl+
'have some default values.'+cl2;

h7=
'L#  is the Longitude of the place / 1 degree, western as negative, default '
 +MyLocation_longitude_string+''+cl+
'D#  the Date (missing values are taken from the system date)'+cl+
'T#  the Time of observation'+cl+
'UT# is the t-UTC difference / 1 h (default 1 from Nov. to March, else 2)'+cl;
h8=
'AL# distance between the shadow tip and the foot of its source'+cl+
'   (ALong the plane, default 1)'+cl+
'P#  how much is the shadow tip nearer to the plane then its source' +cl+
'   (default 0, i.e. Sun lies just within the plane)'+cl;
h9=
'Left/Right - whether the Sun was to the left or right (default Left)'+cl2+

'When the sun shines onto the plane almost tangentially, you can measure'+cl+
'AL and P in such a way: press two matchboxes to the plane some one meter'+cl;
h9a=
'apart and so, that the first bow casts shadow on the second one.'+cl2+

'In another cases, when the shadows are not that long, but still changing'+cl+
'their length with time sufficiently, press a pair of rectangles to the plane'+cl;
h9b=
'so that they stay perpendicular to the plane. One rectangle should have'+cl+
'a circular hole near its end, and should face the Sun.'+cl+
'Mark by the pencil the foot of the shadow-casting rectangle'+cl;

h9c=
'and the centre of spot of light cast through the hole.'+cl+
'Then note the time, and measure the distance (AL) of the spot from the'+cl+
'projection of the hole (it lies inside the outline of foot of rectangle).'+cl;
h9d=
'The quantity P is the distance of the centre of hole from the plane (given'+cl+
'by the rectangle dimensions).'+cl2;

h10=
'M#  9 (default) or 24 needle printer for hardcopy; Ctrl+PrintScreen or P'+cl+
' makes hardcopy to "prn", any # to file named Sundial.e9 (or *.e24)'+cl+
' in the TEMP directory (M1 asks for a CS212-25 printer)'+cl;

h105=
'PT output file with a LaTeX version of the screen (default screen.tex)'+cl+
'PP -     "     -    an Encapsulated PostScript -  "  -     screen.eps)'+cl+
' LA prepend the above eps file with a prologue for a landscape print'+cl;
h11=
'In the resulting plot, the pole is shown as a circle, the cross gives'+cl+
'  the projection of the "gnomon", the shadow-casting point, of the sundial,'+cl+
'  the line on the right side is its distance from the sundial.'+cl;
h12=
'The "hour lines" are given for each degree of the declination of the Sun,'+cl+
'the "hyperboles" are given for solar ecliptical longitude = n * 30 degrees'+cl+
'                           and for each 10 minutes.'+cl;
h13=
'The "8" shows the difference of the true and mean solar time - the arrow'+cl+
'  connects Oct. 11 and 21 (dots are shown for 1.,11. and 21. of each month)'+cl2;

h14=
'E[<name>[:<cf_name>]  shows an exposure table afterwards,'+cl+
'  <name> is a name of file to which the information should be appended,'+cl+
'          <cf_name> is a name of file of the same structure, which should'+cl;
h15=
'          serve as a comparison (default comparison is for'+cl+
'          45 degrees over south at latitude '
 +MyLocation_latitude_string+' degrees).'+cl+
' This is for evaluation of solar collectors and windows. Further parameters'+cl;
h15a=
' for this purpose are:'+cl+
'AB# absorption of radiation in glass at normal incidence (default 0.10)'+cl+
'N#  refractive index of glass (default 1.5, as an average for solar spectrum)'+cl;
h15b=
'Z#[mag|%]  is extinction of solar radiation for Sun in zenith'+cl+
'    (default 0.3 mag, i.e. 24%; default unit is 1 magnitude for #<2'+cl+
'                                and 1 cmag otherwise)'+cl2+

'Q   batch mode: do not wait for a pressed key at the end'+cl;

h16=
'No parameters are obligatory, ? gives this help'+cl2;
h17=
'( (C) Jan Hollan, N.Copernicus Observatory and Planetarium in Brno, 1999-2004;'+cl+
' subject to the GNU General Public License, http://www.gnu.org/copyleft;'+cl+
' source code available at http://astro.sci.muni.cz/pub/hollan/programmes)'+cl;

var
 M,ML,tD_to_Ah: M3;
 V,VA: V3;
 X,Y,AuxReal,AuxReal1,t,Ap,hp : real;
 hour,DecCase,Noon,wrong_pos,jfor: integer;
 IX, IY, PIX, PIY, HIY, SIX: longint;
 index,j : shortint;
 Switch,ch: char;
 Inside, S_In, W_In: boolean;
 outfile: text;
 ZodDec: -3..3;
 text_line:string;
 SumRefnew:array [-3..3] of real;
 SumRel:array [-3..3,0..2] of string[8];
 SumRefYnew:real;
 ro,wR,w2R: real;

const
 rad_large=10;
 ref_file:PathStr='';
 R: real = 1;
 Perp: real = 0;
 Along: real = 1;
 TA: real = 12;
 Loc_Belt: real = 0;
 ShX: integer = 0;
 ShY: integer = 0;
 SmallDesc: boolean = false;
 DescBottom: boolean = false;
 Left: boolean = true;
 Compute_azimuth: boolean = false;
 WFile: boolean=false;
 WFileName : string = 'SunDial.txt';
 North:boolean=false;
 South:boolean=false;  {pole inside the plane}
 YS:byte=94;
 xs:byte=100;
 radius:byte=rad_large;
 Ori_Sub: boolean = true; {origin at the sub-gnomon point}
 AddX: real = 0;
 AddY: real = 0;
 DeclZod:array[-3..3] of real = (-23.4,-20.2,-11.5,0,11.5,20.2,23.4);
 NT: array[-3..3] of byte = (0,0,0,0,0,0,0);
 SumCos: array [-3..3,0..2] of real = ((0,0,0),(0,0,0),(0,0,0),
                                               (0,0,0),
                                       (0,0,0),(0,0,0),(0,0,0));
 Comp_exposure:boolean=false;
 Expo_name:string='';
 n_glass:real=1.5;
 absorb:real=0.1; {absorption in the glass at normal angle of incidence}
 ZenithExt:real=0.3; {mag}
 SumT:word=0;
 SumE:array [0..2] of real=(0,0,0);
 SumRef:array [-3..3] of real =(
{ 1.93,
 2.46,
 3.74,
 4.97,
 5.67,
 5.88,
 5.90);}
1.99,
2.53,
3.85,
5.12,
5.84,
6.06,
6.08);
SumRefY:real= 1647;
               {1599;}
 Sol_Const:real={1.326 much outdated!} 1.366 ;
 BatchMode:boolean=false;
 cylinder:boolean=false;
 shape:string='Planar';
 Help_filename:string='SunDial.hlp';


procedure help_file;
var h_file:text;
begin
 if rt_rewrite(h_file,Help_filename) then
  begin
   write(h_file,h1,cl,
       h2,h2a,h3,h35,h4,h4a,h5,cl,
       h6,h7,h8,h9,h9a,h9b,h9c,h9d,h10,h105,h_svga1,h_svga2,cl,
       h11,h12,h13,h14,h15,h15a,h15b,h16,h17);
   close(h_file);
  end;
{   else help('');}
 halt
end;

procedure help(message:string);
begin
if BatchMode then
 begin
  if (message<>'') and (message<>' ')  then
   begin
    writeln(cl,message,cl);
   end;
  Help_filename:='';
  help_file;
 end
else
 begin   
highvideo;
if message<>'' then
 begin
  writeln(cl,message,cl);
  lowvideo;
   lowvideo;
   write(h0);
   ch:=readkey;
   if UpCase(ch)= 'Q' then halt;
   write(h0clear);
  highvideo;
 end;
writeln(h1);
lowvideo;
write(h2,h2a,h3,h35,h4,h4a,h5,h0);
ch:=readkey;
if UpCase(ch)= 'Q' then halt;
write(h0clear,h6,h7,h8,h0);
ch:=readkey;
if UpCase(ch)= 'Q' then halt;
write(h0clear,h9,h9a,h9b,h9c,h9d,h0);
ch:=readkey;
if UpCase(ch)= 'Q' then halt;
write(h0clear,h10,h105,h_svga1,h_svga2,h0);
ch:=readkey;
if UpCase(ch)= 'Q' then halt;
write(h0clear,h11,h12,h13,h14,h15,h15a,h15b,h0last);
ch:=readkey;
if UpCase(ch)= 'Q' then halt;
write(h0clear,h16,h17);
halt
end;
end;

function Gnom3S(V:V3;R:real;var X,Y:real):boolean;
var aux_boo: boolean; u,sf,cf,f,sigX: real;
begin
 if cylinder and (abs(V[1])>0.001) then
  begin
   u:=sqr(V[3]/V[1]);
   sigX:=V[1]/abs(V[1]);
   if V[3]>0 then  {holds just for ro=R?}
    cf:=(wR + sqrt(w2R - (1+u)*(w2R - u)))/(1+u)
   else
    cf:=(wR - sqrt(w2R - (1+u)*(w2R - u)))/(1+u);
   sf:=sigX*sqrt(1-sqr(cf));
   f:= C2_to_Pf(cf,sf);
   if f>pi then f:=f - 2*pi;
   {f:=arctan(sf/cf);
   if cf<0 then
    if f>0 then f:=pi-f
           else f:=-pi-f;
   }
   X:=ro*f;
   Y:=V[2]*ro*sf/V[1];
   aux_boo:=true;
  end
 else
  aux_boo:=Gnom3(V,R,X,Y);
 if aux_boo and (abs(x)<2E6) and (abs(Y)<2E6) then
         aux_boo:=true else aux_boo:=false;
 Gnom3S:=aux_boo
end;

procedure XYtoIXIY;
begin
        IX:=round(800*(-X))+ShX;
        IY:=round(800*(Y)) +ShY;
end;

procedure PlotPoint;
   begin
    P3S_Vector(-deg_to_rad(t),deg_to_rad(Decl),V);
    VA:=V;
    Mat_Vect(td_to_ah,VA);
    if VA[3]>=0 then
     begin
      Mat_Vect(M,V);
      if Gnom3S(V,R,X,Y) then
       begin
        inside:=true;
        XYtoIXIY;
        if (IX>-1) and (IX<MaxX) and (IY>-1) and (IY<MaxY) then
         {gm_circle_full(IX,IY,radius);}
          gm_star(IX,IY,radius,radius,0,false,'');
       end
      else inside:=false
     end
    else inside:=false
   end;

procedure PlotWPoint;
begin
 PlotPoint;
 if inside then
  begin
   if text_line<>'' then write(outfile,text_line);
   text_line:='';
   write(outfile,-round(X*100+AddX):5,-round(Y*100+AddY):5)
  end
 else if text_line<>'' then text_line:=text_line+'          ';
end;

function Inside_screen:boolean;
begin
 XYtoIXIY;
 Inside_Screen:= (IX>-1) and (IX<=MaxX) and (IY>-1) and (IY<=MaxY)
end;

Procedure SolstIX(SigD,SigT:shortint;var _In:boolean);
begin
    Decl:=deg_to_rad(SigD*23.4);
    V[1]:=-sin(Decl)*sin(fi)/cos(fi);
    V[3]:=sin(Decl);
    V[2]:=SigT*sqrt(1-sqr(V[1])-sqr(V[3]));
    VA:=V;
    Mat_Vect(M,V);
    if Gnom3S(V,R,X,Y) then
     begin
      _In:=true;
      if Inside_screen then
     end
    else
     begin
      _In:=false;
      if V[1]>0 then IX:=0 else IX:=MaxX
     end;
end;

procedure Analema_point;
begin
   Sid_time_loc(SM,JD,11,deg_to_rad(15));
   Ecl_to_Equ(Sun_long(JD),RA,Decl);
   t:=SM*15-RA+Loc_Belt+(TA-12)*15;
   PlotPoint
end;

function A_time(hp,Along,Perp:real; left:boolean; var Ap:real):boolean;
 {Azimuth from hp
    (angular height of the vector pointing normally from the plane to the sky)
 and time, when the shadow cast to a distance of
 >Along<  is  >Perp<  nearer to the plane then its source.
    YYYY,MM,ND and UT, as well as Fi and Lambda
    are to be set prior to calling this procedure.}
var
 C,ApDown,ApHalf,Dif_C,Dif_CD : real;
 NoCHSig: boolean;
begin
 A_time:=true;
 Sid_time_loc(SM,JD,UT,Lambda);
 Ecl_to_Equ(Sun_long(JD),RA,Decl);
 A_h_PA(SM,RA,Decl,A,h,PA,Air_mass);
 A:=deg_to_rad(A);                   {Azimuth of the Sun}
 P3S_Vector(-A,deg_to_rad(h),VA);        {VA is the Sun vector}
 C:= Perp / sqrt(sqr(Along)+sqr(Perp)); {scalar product of the plane's normal
                                     skyward vector and the Sun vector}
 P3S_Vector(-A,hp,V);    {starting from solar azimuth}
 Dif_CD:=scal_prod(VA,V)-C;
 if abs(Dif_CD)<1E-3 then
  begin
   Ap:=A;
   writeln(cl,
'The Sun was almost in the azimuth of the plane. In this situation,',cl,
'the azimuth determination from the length of the shadow is rather inaccurate.',cl,
'Please, make a new measurement in another time, and consider the present value',{);}
{   writeln(}
'to be just a very crude estimate.',cl+cl,

'(press any key to proceed with the program)');
   ch:=readkey;
   exit
  end;

 if left then A:=A-pi;
 P3S_Vector(-A,hp,V);    {... or the opposite one}
 Dif_CD:=scal_prod(VA,V)-C;
 if abs(Dif_CD)<1E-3 then
  begin
   Ap:=A;
   writeln(cl,
'The Sun was almost in the opposite azimuth as that of the plane.',cl,
'In this situation, the azimuth determination from the length of the shadow',cl,
'is rather inaccurate.',cl,
'Please, make a new measurement in another time, and consider the present value',{);}
{   writeln(}
'to be just a very crude estimate.',cl+cl,

'(press any key to proceed with the program)');
   ch:=readkey;
   exit
  end;

 jfor:=0;
 NoChSig:=true;
                                  {finding the interval of intersecion}
 while (jfor<180) and NoChSig do     {of the two (small) circles}
  begin
   inc(jfor);
   Ap:=A+deg_to_rad(jfor);
   P3S_Vector(-Ap,hp,V);
   Dif_C:=scal_prod(VA,V)-C;
   if Dif_C/Dif_CD < 0 then NoCHSIg:=false
  end;
 if NoChSig then begin A_time:=false; exit end;

 ApDown:=A+deg_to_rad(jfor-1);  {the intersection by an iterative halving}
 P3S_Vector(-ApDown,hp,V);
 Dif_CD:=scal_prod(VA,V)-C;
 if Dif_CD<>0 then
 repeat
  ApHalf:=(Ap+ApDown)/2;
  P3S_Vector(-ApHalf,hp,V);
  Dif_C:=scal_prod(VA,V)-C;
  if Dif_C/Dif_CD < 0 then
   Ap:=ApHalf
  else
   begin
    ApDown:=ApHalf;
    Dif_CD:=Dif_C
   end;
 until abs(Dif_C)<1E-5
 else Ap:=ApDown;
 {Ap:=(Ap+ApDown)/2;}
end;

 function Transmitted:real;
  var
   Ai:real; {angle of incidence}
   Ar:real;     {angle of refracted ray}
   Sa:real;
   Da:real;
   S,T: real;
  begin
        Ai:=arctan(sqrt(1-sqr(V[3]))/V[3]);
        Ar:=sin(Ai)/n_glass;
        Ar:=arctan(Ar/sqrt(1-Sqr(Ar)));
         {Now, both Ai and Ar are in radians.}
        Sa:=Ai+Ar;
        Da:=Ai-Ar;
	if Ai<0.01 then {Sa could be zero, 
	                 so Ai should be discarded from the fraction}
	 begin
        S:=sqr((1-1/n_glass)/(1+1/n_glass));
        T:=S;
	 end
	else
	 begin
        S:=sqr(sin(Da)/sin(Sa));
        T:=sqr(sin(Da)*cos(Sa)/(sin(Sa)*cos(Da)));
	 end; 
        Transmitted:=((1-absorb)/cos(Ar))
                     *((1-T)/(1+T)+(1-S)/(1+S))/2;
  end;

procedure Comp_Ex;
begin      
      if inside and Comp_exposure then
       begin
        inc(NT[ZodDec]);
        auxreal:=Sol_Const*V[3]
	 * exp(-L10 * ZenithExt*1.014/(0.014+VA[3]) / 2.5 )
	 ;
        auxreal1:=transmitted;
        SumCos[ZodDec,0]:=SumCos[ZodDec,0]+auxreal*auxreal1;
        SumCos[ZodDec,1]:=SumCos[ZodDec,1]+auxreal;
                  {not the reflected part}
        SumCos[ZodDec,2]:=SumCos[ZodDec,2]+auxreal*sqr(auxreal1);
      end; 
end;


begin     {of the main program SunDial}
FileMode:=0; {will open the files by reset as Read Only
               -- they can be then marked so even if not on a CD;
               here it enables to have *.bgi and *.chr read-only}
 Main_program:='SunDial';
 Ap:= 0;
 hp:= 0;
 Decl:=0;
 if paramcount > 0 then  for jfor:=1 to paramcount do
  begin
   wrong_pos:=0;
   Par:=paramstr(jfor);
   if (Par[1]='/') or (Par[1]='-') then Par:=copy(Par,2,length(par)-1);
   Switch:=UpCase(Par[1]);
   if length(Par)=1 then
    case Switch of
     'S': SmallDesc:=true;
     'B': begin SmallDesc:=true; DescBottom:=true end;
     'O': ori_sub:=false;
     '?','-': Help(''); {- added to respect unix --help even without the help}
    end;
   if (length(Par)>1) or not (Switch in ['S','B','O']) then
     begin
      Par:=copy(Par,2,length(Par)-1);
      val(Par,AuxReal,wrong_pos);
      case Switch of
       'A': if length(par)>1 then
             case UpCase(Par[1]) of
             'L':
               begin    {distance of the shadow from its source along the plane}
                Par:=copy(Par,2,length(Par)-1);
                val(Par,Along,wrong_pos)
               end;
             'B':if length(Par)>1 then
               begin
                absorb:=par2r(par);
               wrong_pos:=0;
               end;
             else
              Ap:= deg_to_rad(AuxReal-180);
             end
            else
             Ap:= deg_to_rad(AuxReal-180);
       'B': if length(par)>0 then
             begin
              case UpCase(Par[1]) of
               'C','R':
                    if length(par)>1 then
                     if Par[2]='-' then Draw_picture_border:=' '
                     else Draw_picture_border:=Par[1];
                    else Draw_picture_border:=Par[1];
              end;
              wrong_pos:=0;
             end;
       'C': if wrong_pos=0 then
             begin
              ro:=AuxReal;
              cylinder:=true;
             end;
       'D': if Par<>'' then
             if UpCase(Par[1])='L'then
              begin        {difference of the desired time from the local one}
               Par:=copy(Par,2,length(Par)-1);
               val(Par,AuxReal,wrong_pos);
               Loc_Belt:=AuxReal/4
              end
             else
              begin
               Par:='D'+Par;
               if not par_date then help(' ');
               wrong_pos:=0;
              end;
       'E': begin
             Comp_exposure:=true;
             if Par<>'' then
              begin
               Expo_name:=ItemStrD(':',1,Par);
               if ItemCountD(':',Par)>1 then
                ref_file:=ItemStrD(':',2,Par);
              end;
             wrong_pos:=0;
            end;
       'F': begin
             Par:='F'+Par;
             if not Par_Fi then help(' ')
            end;
       'H': if Par='' then help_file
            else hp:= deg_to_rad(AuxReal);
       'L': if wrong_pos>0 then
             begin       {Sun was to the left}
              if length(Par)>0 then
               case  UpCase(Par[1]) of
               'A': LandscapeA4:=true;  {added 2002-08-01, more paper sizes
                                     could be included}
               else Left:=true
               end;
               wrong_pos:=0
             end
            else
             Lambda:=deg_to_rad(AuxReal);
       'M': mode:=round(AuxReal);
       'N': if length(Par)>0 then
             begin auxreal:=ss2r(Par,1); if auxreal>1 then n_glass:=auxreal;
             end;
       'P': if Par<>'' then
             case UpCase(Par[1]) of
              'T':
                begin
                 GM_mode:=2;
                 GM_OuFile:= copy (Par,2,length(Par)-1);
                 wrong_pos:=0;
                end;
               'P':
                begin
                 GM_mode:=3;
                 GM_OuFile:= copy (Par,2,length(Par)-1);
                 wrong_pos:=0;
                end;
             else
              Perp:=AuxReal;
             end;
       'Q': begin BatchMode:=true; wrong_pos:=0 end;
       'R': if wrong_pos>0 then
             begin       {Sun was to the right}
              Left:=false;
              wrong_pos:=0
             end
            else
             R:=AuxReal;
       'S': if length(par)>2 then
             if __cvtstr(copy(Par,1,3),_to_upcase_str)='VGA' then
              begin
               Par_SVGA;
               wrong_pos:=0;
              end;
       'T': if par<>'' then
             if UpCase(Par[1])='A'then       {analema time}
              begin
               Par:=copy(Par,2,length(Par)-1);
               val(Par,TA,wrong_pos)
              end
             else
              begin
               Par:='T'+par;
               wrong_pos:=0;
               if not par_time then help(' ');
               Compute_azimuth:=true
              end;
       'U': begin
             Par:='U'+Par;
             if not Par_UT then help(' ');
             wrong_pos:=0;
            end;
       'W': begin WFile:=true; if Par<>'' then WFileName:=Par; Wrong_pos:=0 end;
       'X': if Par<>'' then
             if UpCase(Par[1])='S' then
              val(copy (Par,2,length(Par)-1),XS,wrong_pos)
             else
              ShX:=round(AuxReal);
       'Y': if Par<>'' then
             if UpCase(Par[1])='S' then
              val(copy (Par,2,length(Par)-1),YS,wrong_pos)
             else
              ShY:=round(AuxReal);
       'Z': begin
             if wrong_pos>0 then
              begin
	        AuxReal:=R_S(copy(Par,1,wrong_pos-1));
               case par[wrong_pos] of
	        '%' : AuxReal:=-2.5*ln(1-AuxReal/100)/l10;
	        'c' : AuxReal:=AuxReal/100;
	        'd' : AuxReal:=AuxReal/10;
               end;
              end
	     else if AuxReal>2 then AuxReal:=AuxReal/100; 
             ZenithExt:=AuxReal;
            end;
       else help('"'+Switch+'" is no parameter from the allowed set.'+
               ' See the following instructions:');
      end {of case}
    end;
   if wrong_pos>0 then
     help('Sorry, the parameter "'+Paramstr(jfor)+'" gives no number - '+
             SI(wrong_pos,2)+'-th character is bad.');
  end; {of for}

 if Wfile then if not rt_rewrite(outfile,WFileName) then Wfile:=false;

 if cylinder then
  begin
   wR:=1-1/ro;
   w2R:=sqr(wR);
   ro:=R*ro;
   shape:='Cylindric';
   Comp_exposure:=false;
  end;

 if Compute_Azimuth then
  begin
   UT:= HH+(int(MI)/60)+(int(SS)/3600) - t_UT;
   ND:=DD;
   if UT<0  then dec(ND);
   if UT>24 then inc(ND);
   UT:=Bas_int(UT,24);
   if not A_time(hp,Along,Perp,left,Ap) then
     help('The Sun could not shine onto the plane that way at that moment.');
  end;

 Mat(M, 2, pi/2 - Fi); {Z to zenith}
 tD_to_Ah:=M;
 Mat(ML, 3, -Ap);
 Matl_mat(ML,M);       {X to Ap}
 Mat(ML, 2, pi/2 - hp);
 Matl_mat(ML,M);       {Z to hp}
 Mat(ML, 3, pi/2);
 Matl_mat(ML,M);       {X to be horizontal}
 gm_init(xs,ys);

 gm_comm('Description of the task:');
 if SmallDesc then
  begin
   if DescBottom then IY:=MaxY-TextH-2 else IY:=TextH;
   gm_outxy(MaxX div 2,IY,'A='+SRe(4,1,Bas_int(rad_to_deg(Ap)+180,360))+
           ', h='+SRe(4,1,rad_to_deg(hp))+
           ', Fi='+SRe(4,1,rad_to_deg(Fi))+' degrees'+
           ', Loc.time '+SRe(4,1,Loc_Belt*4)+' min ahead of the shown one.',
           {' (and R='+SRe(4,2,R)+',X='+SI(2,ShX)+',Y='+SI(2,ShY)+
                 ', T='+SRe(4,2,TA)+')'}
           'c',0);

  end
 else
  begin
   IY:=50;
   gm_outxy(32,IY,shape+' sundial for these parameters:','l',0);
   inc(IY,TextH);
   if cylinder then 
    begin
     gm_outxy(32,IY,
     'Radius of the cylinder is'+SRe(7,2,ro/R)+
     ' times the distance of the shadow-casting point from the plane','l',0);
     inc(IY,TextH);
     gm_outxy(32,IY,'The line perpendicular to its tangent plane points at:','l',0)
    end
   else    
    gm_outxy(32,IY,'The line perpendicular to its plane points at:','l',0);
   inc(IY,TextH);
   gm_outxy(32,IY,'   azimuth  '+SRe(4,1,Bas_int(rad_to_deg(Ap)+180,360))+
              ' degrees,','l',0);
   inc(IY,TextH);
   gm_outxy(32,IY,'   angular height  '+SRe(4,1,rad_to_deg(hp))+
              ' degrees,','l',0);
   inc(IY,TextH);
   gm_outxy(32,IY,'The latitude of the place is '+SRe(4,1,rad_to_deg(Fi))+
                   ' degrees','l',0);
   inc(IY,TextH);
   gm_outxy(32,IY,'The local time is '+SRe(4,1,Loc_Belt*4)+
                   ' min ahead of the shown time','l',0);
   inc(IY,TextH);
   gm_outxy(32,IY,'(SunDial /? would give you help.)','l',0);
  end;

 if Wfile then
  begin
   writeln(outfile,shape+' sundial for these parameters:');
   if cylinder then 
    begin
     writeln(outfile,
      'Radius of the cylinder is'+SRe(5,2,ro/R)+
      ' times the distance R of the shadow-casting'+cl+
      '  point from the plane, i.e. radius = ',SRe(7,2,Ro*100));
     writeln(outfile,'The line perpendicular to its tangent plane points at:','l',0)
    end
   else    
    writeln(outfile,'The line perpendicular to its plane points at:');
   writeln(outfile,'   azimuth  '+SRe(4,1,rad_to_deg(Ap)+180)+' degrees,');
   writeln(outfile,'   angular height  '+SRe(4,1,rad_to_deg(hp))+' degrees,');
   writeln(outfile,'The latitude of the place is '+SRe(4,1,rad_to_deg(Fi))+
                   ' degrees');
   writeln(outfile,'The local time is '+SRe(4,1,Loc_Belt*4)+
                   ' min ahead of the shown time');
   writeln(outfile,'The shadow-casting point is at distance R='+SRe(7,2,R*100),
                   ' from the plane',cl);
   if Ori_sub then
    writeln(outfile,'The coordinates origin at the sub-gnomon point,')
   else
    writeln(outfile,'The coordinates origin at the pole,');
   writeln(outfile,' X goes horizontally to the right, Y upwards.',cl);
  end;

 ShX:=(MaxX div 2) + round((MaxX+1)*ShX/100); {to the screen coordinates}
 ShY:=(MaxY div 2) - round((MaxY+1)*ShY/100);

 gm_comm(cl+ '% projection of the centre of the  gnomonic projection');
 gm_comm('and of its distance from the plane as line segments at +-x, +-y:');
 for i:=1 to 2 do V[i]:=0;
                V[3]:=1;
 if Gnom3S(V,R,X,Y) then;
   if Inside_screen then
    gm_line(IX-80,IY,IX+80,IY);
   auxReal:=X;
   X:=X+R;
   if Inside_screen then
    begin
     gm_line(IX,IY,IX+80,IY);
     gm_line(IX,IY-32,IX,IY+32)
    end;
   if cylinder then 
    begin           
     X:=auxReal + ro;
     if Inside_screen then
      begin
       gm_comm(cl+ '% cylinder radius plotted right from the cross:');
       gm_line(IX,IY,IX+80,IY);
       gm_line(IX,IY-32,IX,IY+32)
      end; 
     X:=auxReal - ro;
     if Inside_screen then
      begin
       gm_comm(cl+ '% cylinder radius plotted left from the cross:');
       gm_line(IX-80,IY,IX,IY);
       gm_line(IX,IY-32,IX,IY+32)
      end;
     X:=auxReal+R; 
     gm_comm(cl+cl+'% centre of projection distance from the plane continuation:');
    end;
   X:=X-2*R;
   if Inside_screen then
    begin
     gm_line(IX-80,IY,IX,IY);
     gm_line(IX,IY-32,IX,IY+32)
    end;
   X:=X+R;
   if Inside_screen then
    gm_line(IX,IY-80,IX,IY+80);
   Y:=Y-R;
   if Inside_screen then
    begin
     gm_line(IX,IY,IX,IY+80);
     gm_line(IX-32,IY,IX+32,IY)
    end;
   Y:=Y+2*R;
   if Inside_screen then
    begin
     gm_line(IX,IY-80,IX,IY);
     gm_line(IX-32,IY,IX+32,IY)
    end;

 VA:=V;
 Mat_Vect(M,V);
 if Gnom3S(V,R,X,Y) then
  begin
   if Inside_screen then
    begin
     gm_comm(cl+'% the north pole:');
     gm_circle_empty(IX,IY,48);
    end;
   PIX:=IX; PIY:=IY;
   North:=true;
   if Wfile then
    if Ori_sub then
     writeln(outfile,'Pole X=',-round(X*100):5,' Y=',-round(Y*100):5,cl)
    else
     begin
      writeln(outfile,'Sub Gnomon X=',round(X*100):5,' Y=',round(Y*100):5,cl);
      AddX:=-X*100; AddY:=-Y*100;
     end;
  end;
 V:=VA; V[3]:=-1;
 Mat_Vect(M,V);
 if Gnom3S(V,R,X,Y) then
  begin
   if Inside_screen then
    begin
     gm_comm(cl+'% the south pole:');
     gm_circle_empty(IX,IY,48);
    end;
   PIX:=IX; PIY:=IY;
   South:=true;
   if Wfile then
    if Ori_sub then
     writeln(outfile,'Pole X=',-round(X*100):5,' Y=',-round(Y*100):5,cl)
    else
     begin
      writeln(outfile,'Sub Gnomon X=',round(X*100):5,' Y=',round(Y*100):5,cl);
      AddX:=-X*100; AddY:=-Y*100;
     end;
  end;

 if Wfile then
  begin
   writeln(outfile,'Length of the shaft projection is',round(sqrt(X*X+Y*Y)*100):5,cl);
   writeln(outfile,'Length of the shaft is',round(sqrt(X*X+Y*Y+R*R)*100):5,cl);
  end;
 if (hp < pi/2 -0.1) and not cylinder then
  begin
    gm_comm(cl+'% horizontal plane:');
    V[1]:= 0;
    V[2]:=-sin(hp);
    V[3]:= cos(hp);
    if Gnom3S(V,R,X,Y) then;  {of course it is true}
    if Inside_screen then;
    if Wfile then
     writeln(outfile,'The h=0 line has Y=',-round(Y*100+AddY):5,cl);
    HIY:=IY;

    SolstIX(1,1,S_In);
    SIX:=IX;
    SolstIX(-1,1,W_In);
    if (not(SIX=IX)) and (S_In or W_In) then
     gm_line(SIX,HIY,IX,HIY);
    SolstIX(1,-1,S_In);
    SIX:=IX;
    SolstIX(-1,-1,W_In);
    if (not(SIX=IX)) and (S_In or W_In) then
     gm_line(SIX,HIY,IX,HIY);
  end;

 YYYY:=2000;
 Radius:=9;
 if Wfile then writeln(outfile,
  'Analema for the plotted time T='+SRe(4,2,TA)+':'+cl+'MM-DD   X    Y');
 gm_comm(cl+'% analema:');
 for MM:=1 to 12 do for jfor:=0 to 2 do
  begin
   ND:=10*jfor+1;
   Analema_point;
   if WFile then
    begin
     writeln(outfile,MM:2,'-',ND:2,-round(X*100+AddX):5,-round(Y*100+AddY):5);
     if jfor=2 then writeln(outfile);
    end;
  end;
 if Wfile then writeln(outfile);

 MM:=10;
 Radius:=1;
 gm_comm(cl+'% arrow to Oct.21:');
 for jfor:=1 to 2  do
 begin
   ND:=1+10*jfor;
   Analema_point;
   if Jfor=1 then begin SIX:=IX; HIY:=IY end;
 end;
 gm_line(SIX,HIY,IX,IY);
 SIX:=IX; HIY:=IY;
 ND:=19;
 Analema_point;
 t:=t+0.3;
 PlotPoint;
 gm_line(SIX,HIY,IX,IY);
 t:=t-0.6;
 PlotPoint;
 gm_line(SIX,HIY,IX,IY);

 Radius:=rad_large;
 gm_comm(cl+'% points (-11, -10,... to 12h; -23, -22,... to 23 degrees):');
 for hour:=-11 to 12 do
  begin
  gm_comm(cl+'% hour'+SI(4,hour+12)+':');
  t:=hour*15+Loc_Belt;
  for DecCase:=-23 to 23 do {lines}
   begin
    Decl:=DecCase;
    PlotPoint;
   end;
  if (North or South) and not cylinder then
   begin
    if North then Decl:=23.4 else Decl:=-23.4;
    P3S_Vector(-deg_to_rad(t),deg_to_rad(Decl),V);
    VA:=V;
    Mat_Vect(td_to_ah,VA);
    if VA[3]>=0 then   {seems to mean: Sun in a positive angular height}
     begin
      Mat_Vect(M,V);
      if Gnom3S(V,R,X,Y) then
       begin
        inside:=true;
        XYtoIXIY
       end
      else inside:=false
     end
    else inside:=false;
    if not inside then
     begin
      P3S_Vector(-deg_to_rad(t),arctan(-cos(deg_to_rad(t))*cos(Fi)/sin(Fi)),V);
      Mat_Vect(M,V);
      if Gnom3S(V,R,X,Y) then
       begin
        inside:=true;
        XYtoIXIY
       end
      else inside:=false
     end;
    if inside then
     begin
      gm_comm('line from pole:');
      gm_line(IX,IY,PIX,PIY)
     end
   end;
  for DecCase:=-1 to 1 do if DecCase<>0 then {description}
   begin
    if DecCase=1 then Noon:=13 else Noon:=12;
    if North then if DecCase=-1 then Decl:={-25}-33 else Decl:=33
    else
     if South then if DecCase=-1 then Decl:=-33 else Decl:={25}33
     else Decl:=DecCase*{26}33;
    P3S_Vector(-deg_to_rad(t),deg_to_rad(Decl),V);
    VA:=V;
    Mat_Vect(td_to_ah,VA);
    if VA[3]>-0.3 {description can be over horizontal line} then
     begin
      Mat_Vect(M,V);
      if Gnom3S(V,R,X,Y) then
        if Inside_screen then
         begin
          gm_comm('description:');
          if hour+Noon>9 then W:=2 else W:=1;
          gm_outxy(IX,IY {+ texth div 3},SI(W,hour+Noon),'w',0)
         end
     end
   end
  end;

 gm_comm(cl+
'% hyperboles for declinations 23.4, 20.2, 11.5, 0, -11.5, -20.2, -23.4 degrees:');
 if WFile then
 begin
  writeln(outfile,
'Decl.  -23.4     -20.2     -11.5      0         11.5      20.2      23.4');
  write(outfile,
'hh:mm');
  for hour:=1 to 7 do write(outfile,      {hour is just column No. here}
     '  X    Y  ');
  writeln(outfile);
   for hour:=-71 to 72 do {each 10 minutes}
   begin
    t:=hour*2.5+Loc_Belt;
    text_line:=SI(2,(hour+72) div 6)+':'+SI(2,((hour+72) mod 6)*10);
    gm_comm(text_line);
    for ZodDec:=-3 to 3 do
     begin
      Decl:=DeclZod[ZodDec];
      PlotWPoint;
      Comp_Ex;
     end;
    if text_line='' then
     begin
      writeln(outfile);
      if (hour+72) mod 6 = 5 then writeln(outfile);
     end;
   end
 end
 else
   for hour:=-71 to 72 do {each 10 minutes}
   begin
    t:=hour*2.5+Loc_Belt;
    gm_comm(SI(2,(hour+72) div 6)+':'+SI(2,((hour+72) mod 6)*10));
    for ZodDec:=-3 to 3 do
     begin
      Decl:=DeclZod[ZodDec];
      PlotPoint;
      Comp_Ex;
     end;
   end;

 if not BatchMode then
  begin
   ch:=readkey;
   if ch=#0 then
    begin
     ch:=readkey;          {PrintScreen}
     if ch=#114 then PriGScr('',0,mode)
    end
   else
    if UpCase(ch)='P' then PriGScr('',0,mode)
    else
     begin
      if ord(ch) in [ord('0')..ord('9')] then
       PriGScr('SunDial',ord(ch)-ord('0'),mode);
     end;
    ch:=' ';
  end;
 if UpCase(ch)<>'G' then gm_close;   {G leaves the map on the screen}
 if WFile then close(outfile);
 if Paramcount=0 then help('');
 if Comp_exposure then
  begin
   if ref_file<>'' then
    Wfilename:=fsearch(ref_file,'')
   else
    Wfilename:='';
   index:=0;
   if Wfilename<>'' then       {taking reference values from ref_file:}
    begin
     if rt_reset('',outfile,ref_file) then index:=-3;
     repeat
      readln(outfile,text_line);
      if itemcount(text_line)>3 then
       begin
        val(ItemStr(1,text_line),auxreal,wrong_pos);
        if (wrong_pos=0) or (index=4) then
         begin
          val(ItemStr(2,text_line),auxreal,wrong_pos);
          if wrong_pos=0 then
           begin
            val(ItemStr(3,text_line),auxreal,wrong_pos);
             if wrong_pos=0 then
              if index=4 then
               begin
                SumRefYnew:=auxreal;
                inc(index);
               end
              else
               begin
                SumRefnew[index]:=auxreal;
                inc(index);
               end
           end
         end
       end;
     until eof(outfile) or (index=5);
     close(outfile);
     if index=5 then  {setting new reference values:}
      begin
       for ZodDec:=-3 to 3 do SumRef[ZodDec]:=SumRefnew[ZodDec];
       SumRefY:=SumRefYnew;
      end
    end;
   for ZodDec:=-3 to 3 do    {relative values as strings; sums}
    begin
     for j:=0 to 2 do
      if SumRef[ZodDec]>0 then
       SumRel[ZodDec,j]:=SRe(8,2,SumCos[ZodDec,j]/(6*SumRef[ZodDec]))
      else
       SumRel[ZodDec,j]:='infinity';
     if abs(ZodDec)=3 then
      begin
       SumT:=SumT+NT[ZodDec]*5;
       for j:=0 to 2 do
        SumE[j]:=SumE[j]+SumCos[ZodDec,j]*5;
      end
     else
      begin
       SumT:=SumT+NT[ZodDec]*10;
       for j:=0 to 2 do
        SumE[j]:=SumE[j]+SumCos[ZodDec,j]*10;
      end
    end;
   if Expo_name<>'' then
    begin
     if fsearch(Expo_name,'')<>'' then
      Wfilename:=cl+'    -------------------------------'+cl2
     else
      Wfilename:='';
    end;

   if rt_append('',outfile,Expo_name) then;

   writeln(outfile,
'Daily solar exposure of a plane sized one square meter for these parameters:',cl,
' The normal to the plane points to',cl,
'   azimuth  '+SRe(5,1,rad_to_deg(Ap)+180)+' degrees ',
       ' (South is at 180 degrees, West at 270 or -90),',cl,
'   height   '+SRe(5,1,rad_to_deg(hp))+' degrees ',
       ' (Zenith is at 90 degrees),',cl,
' latitude of the place is '+SRe(4,1,rad_to_deg(Fi))+' degrees,',cl,
' sky is CLEAR (all the time!), sunshine diminished',
 SRe(3,0,(1-exp(-0.4*L10*ZenithExt))*100),' %  per unit air mass,',cl,
' glass properties: n =',n_glass:4:2,', absorption = ',absorb:4:2,' (at normal incidence).');

   if index=5 then
    writeln(outfile,cl,
'Reference exposures for column "relative" were taken from file: ',cl,
 ref_file,cl)
   else
    writeln(outfile,
'A reference plane points to 45 degrees over south',cl,
' at latitude of '
 +MyLocation_latitude_string+' degrees (positive=northern), is single-glazed'
,cl,
' and its insolation is diminished by 24 % per unit air mass.',cl
);

   writeln(outfile,
'Declination  Insol-   single glazing       no glazing       double glazing',cl,
'of the Sun   ation  exposure    relat.    Expos.  relat.    Expos.  relat.' ,cl,
'/ 1 degree   time    / 1 kWh             / 1 kWh           / 1 kWh',cl,
'             /1 h');

   for ZodDec:=-3 to 3 do
    begin
     write(Outfile,
DeclZod[ZodDec]:5:1,NT[ZodDec]/6:12:1,' ');
     for j:=0 to 2 do write(OutFile,
       SumCos[ZodDec,j]/6:10:2,SumRel[ZodDec,j]);
     writeln(OutFile);
    end;

   write(outfile,
'Yearly_Sum:',SumT:6,' ');
   for j:=0 to 2 do write(OutFile,
SumE[j]:10:0,SumE[j]/SumRefY:8:2);
writeln(outfile); 
{end of line and an additional blank line
 two lines are needed to prevent doubling the last 
 line from an output of sundialt.php} 
writeln(outfile);

   close(outfile);

  end; {of if Comp_exposure}
end.