{~ Angular and matrix operations, Opening text files ~}

unit angles_o;          {various utilities useful to me, Jan Hollan:
            Opening text files reliably,
     Angular operations, including 3d-rotations }
(*  Copyright (C) 1999 Jan Hollan;
  by "program" this unit is further meant.

    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.
*)
interface
uses dos,
     str_num;
type c3d = 1..3;
     M3 = array [1..3,1..3] of real;
     V3 = array [1..3] of real;
var coord,i,j,k,W: C3d;
 Ma_Ah_to_td:M3;
const
 V_unit_base:array[1..3] of V3 = ((1,0,0),(0,1,0),(0,0,1));

function rt_reset(Stand_Loc:string;var f1: text;f:PathStr):boolean;



      {reliable text assign & reset, false if unsuccessful}
function rt_reset_suf(Stand_Loc:string;var f1: text;f:PathStr;suf:ExtStr):boolean;
      {reliable text assign & reset, false if unsuccessful.
       Like rt_reset, but tries suffix suf if unsuccessful with f alone.}
function rt_rewrite(var f1: text; f:string):boolean;
      {reliable text assign & rewrite, false if unsuccesful}
function rt_append(Stand_Loc:string;var f1: text;f:PathStr):boolean;
      {reliable text assign & append or rewrite, false if unsuccesful}

FUNCTION deg_to_rad(A:REAL):REAL; {degrees to radians}
FUNCTION rad_to_deg(A:REAL):REAL;  {radians to degrees}

Function Bas_int(A:REAL;P:integer):real; {adding n*P to A to reach <0,P)}

function scal_prod(Vf,Vs:V3):real;
procedure vect_prod(Vf,Vs:V3;var Vr:V3);
              {angles in radians, when not explicitly said that in degrees}

procedure vect_prod_dir(l1,b1,l2,b2:real; var lp,bp:real);
              {just direction of vector product}
procedure vect_prod_dirD(l1,b1,l2,b2:real; var lp,bp:real);
              {just direction of vector product, all angles in degrees}
function  C2_to_Pf (x,y:real):real;  {cartesian 2-d system to angle}
function  C2_to_Pdf(x,y:real):real;  {cartesian 2-d system to angle in degrees}
procedure C2_to_Pd (x,y:real; var r,p:real);
                                {cartesian 2-d system to polar one (degrees)}
procedure P2_to_C(r,p:real; var x,y:real);
                              {polar 2-d system           to cartesian one}
procedure P2d_to_C(r,p:real; var x,y:real);
                              {polar 2-d system (degrees) to cartesian one}
procedure C3_to_P(x,y,z:real; var r,l,f:real);
                              {cartesian 3-d system to polar one (radians)}
procedure C3_to_Pd(x,y,z:real; var r,l,f:real);
                                                                {(degrees)}
procedure P3_to_C(r,l,f:real; var x,y,z:real);
                              {polar 3-d system to cartesian one}
procedure P3d_to_C(r,l,f:real; var x,y,z:real);
                                                                {(degrees)}

procedure Vector_P3s(V:V3; var l,f:real);
                              {3-vector to polar sky}
procedure P3s_Vector(L,F:real; {like RA and Decl}
                 var V:V3 {V[3] to max F, V[1] to 0 L});

procedure mat(var M:M3 {matrix}; coord:c3d; omega:real{in radians});
          {transformation matrix for Rotation by Omega around Coord}
procedure mat_vect(M:M3; var V:V3 {vector});
          {vector:= matrix * vector}
procedure matl_mat(ML:M3; var M:M3);
          {matrix:= matrix_left * matrix}
procedure matV(var M:M3 {matrix}; VNP_3,V0L_1:V3);
          {transformation matrix into new base}
procedure matVinv(var M:M3 {matrix}; VNP_3,V0L_1:V3);
          {transformation matrix from the new base to the standard one}
function  Gnom3 (V:V3 {vector};
                 R:real {distance to the projection plane (along V[3])};
                 var X,Y: real {gnomonic coordinates} )
               : boolean {true if V[3]>0};

procedure Ah_to_td(Pole_rot:real;var L,F:real);
                      {lefthand systems like horizontal->equatorial
                    (then Pole_rot=-(pi/2-Fi), backwards pi/2-Fi), radians}
procedure Ah2td(var L,F:real); {like ah_to_td, but uses its Matrix ...}
procedure ra_d2l_b(Pole_rot:real;var L,F:real);
                  {righthand systems equatorial->ecliptical-like, radians,
                   pole_rot = epsi_rad for really eq.->ecl. transformation,
                             -epsi_rad for        ecl.->equatorial transf.}
function object_phase_angleD(l_helioc,b_helioc,l_geoc,b_geoc:real):real;
  {ecliptical coordinates of the object from the Sun and the Earth
   output and input in degrees, ouput is angle Sun-object-Earth}

procedure ang_dist_posi_angD2(ra1_d,de1_d,ra0_d,de0_d,ra2_d,de2_d:real;var pa,ad: real);
 {angular distance and position angle of the second point with respect
  to the first one; position angle is measured counterclockwise from
  the colur going through ra0,de0 (in astronomy usually the current Earth
  celestial north pole with de0_d=90;
  righthand system like equatorial coordinates, output and input in degrees}

procedure ang_dist(ra1_d,de1_d,ra2_d,de2_d:real;var ad: real);


implementation

const li=2147483647;

function rt_reset(Stand_Loc:string;var f1: text;f:PathStr):boolean;
 {if no directory is given
  in f and that file is not found in the current one, searching
  through the StandLoc (DirString: series of directories delimited by ;).
  If still not found writing a message on the screen.}
var FoundName: PathStr; FDir: DirStr; FName: NameStr; FExt: ExtStr;
begin
 rt_reset:=false;
 FSplit(f,FDir,FName,FExt);
 if FDir<>'' then
  FoundName:=FSearch(f,'')
 else
  begin
   FoundName:=FSearch(f,'');
   if FoundName='' then
    FoundName:=FSearch(f,Stand_Loc)
  end;
 if FoundName<>'' then
  begin
   assign(f1,foundname);
   reset(f1);
   rt_reset:=true;
  end
 else
  begin
    writeln;
    write('There is no file >>');
    write(f);
    writeln('<<. Spell the name correctly, please.');
    rt_reset:=false;
  end
end;

function rt_reset_suf(Stand_Loc:string;var f1: text;f:PathStr;suf:ExtStr):boolean;
 {if no directory is given
  in f and that file is not found in the current one, searching
  through the StandLoc (DirString: series of directories delimited by ;).
  Tries suffix >suf< if unsuccessful with >f< alone.
  If still not found writing a message on the screen.}
var FoundName: PathStr; FDir: DirStr; FName: NameStr; FExt: ExtStr;
begin
 rt_reset_suf:=false;
 FSplit(f,FDir,FName,FExt);
 if suf[1]<>'.' then suf:='.'+suf;
 if FDir<>'' then
  begin
   FoundName:=FSearch(f,'');
   if (FoundName='')
   and (length(FExt)<2) then
    FoundName:=FSearch(FDir+FName+suf,'');
  end
 else
  begin
   FoundName:=FSearch(f,'');
   if FoundName='' then
    FoundName:=FSearch(f,Stand_Loc);
   if (FoundName='')
   and (length(FExt)<2) then
    FoundName:=FSearch(FName+suf,Stand_loc);
  end;
 if FoundName<>'' then
  begin
   assign(f1,foundname);
   reset(f1);
   rt_reset_suf:=true;
  end
 else
  if FExt<>'' then  {to be able to try a default suffix next time}
   begin
    writeln;
    write('There is no file >>',f,'<<');
    if (length(FExt)<2) then write ('nor >>',FDir+FName+suf,'<<.');
    writeln(' Spell the name correctly, please.');
    rt_reset_suf:=false;
   end
end;

function rt_rewrite(var f1: text; f:string):boolean;
begin
 {$I-}
 assign(f1,f);
 rewrite(f1);
 if ioresult<>0 then
  begin
   writeln;
   write('The output file >>');
   write(f);
   writeln('<< cannot be opened. Wrong name (or disk full?).');
   rt_rewrite:=false
  end
 else rt_rewrite:=true  ;
 {$I+}
end;

function rt_append(Stand_Loc:string;var f1: text;f:PathStr):boolean;
 {if no directory is given
  in f and that file is not found in the current one, searching
  through the StandLoc (DirString: series of directories delimited by ;).
  If still not found writing a message on the screen.}
var FoundName: PathStr; FDir: DirStr; FName: NameStr; FExt: ExtStr;
begin
 rt_append:=false;
 FSplit(f,FDir,FName,FExt);
 if FDir<>'' then
  FoundName:=FSearch(f,'')
 else
  begin
   FoundName:=FSearch(f,'');
   if FoundName='' then
    FoundName:=FSearch(f,Stand_Loc)
  end;
 if FoundName<>'' then
  begin
   assign(f1,foundname);
   append(f1);
   rt_append:=true;
  end
 else
  begin
   {$I-}
   assign(f1,f);
   rewrite(f1);
   if ioresult<>0 then
    begin
     writeln;
     write('The output file >>');
     write(f);
     writeln('<< cannot be opened. Wrong name (or disk full?).');
    end
   else rt_append:=true  ;
   {$I+}
  end
end;

FUNCTION deg_to_rad(A:REAL):REAL; {degrees to radians}
 BEGIN
  deg_to_rad:=PI*A/180
 END;

FUNCTION rad_to_deg(A:REAL):REAL;  {radians to degrees}
 BEGIN
  rad_to_deg:=180*A/PI
 END;

Function Bas_int(A:REAL;p:integer):real; {adding n*P to reach <0,P)}
 BEGIN
  if not ( (A>-li) and (A<li) ) then
   begin
    writeln('an argument of Bas_Int was outside ',-li,' to ',li,' limit: ',A);
    {Bas_Int:=0;} halt;
   end
  else
   begin
    A:=frac(A) + (trunc(A) mod p);
    if A<0 then
    A:=A+p;
    Bas_int:=A;
   end;
 END;

function scal_prod(Vf,Vs:V3):real;
var AuxReal:real;
begin
 AuxReal:=0;
 for i:=1 to 3 do AuxReal:=AuxReal+Vf[i]*Vs[i];
 scal_prod:=AuxReal
end;

procedure vect_prod(Vf,Vs:V3;var Vr:V3);
begin
 Vr[3]:=Vf[1]*Vs[2]-Vf[2]*Vs[1];
 Vr[2]:=Vf[3]*Vs[1]-Vf[1]*Vs[3];
 Vr[1]:=Vf[2]*Vs[3]-Vf[3]*Vs[2];
end;

function C2_to_Pf(x,y:real):real;
var p:real;
begin
 if x=0 then if y>=0 then p:=pi/2 else p:=3*pi/2
 else
  begin
   p:=arctan(y/x);
   if  x<0  then p:=p+pi
   else if y<0 then p:=p+2*pi
  end;
 C2_to_Pf:=p;
end;

function C2_to_Pdf(x,y:real):real;
begin
 C2_to_Pdf:=rad_to_deg(C2_to_Pf(x,y));
end;

procedure C2_to_Pd(x,y:real; var r,p:real);
begin
 r:=sqrt(sqr(x)+sqr(y));
 p:=C2_to_Pdf(x,y);
end;

procedure P2_to_C(r,p:real; var x,y:real);
                              {polar 2-d system           to cartesian one}
begin
 x:=r*cos(p);
 y:=r*sin(p);
end;


procedure P2d_to_C(r,p:real; var x,y:real);
begin
 P2_to_C(r,deg_to_rad(p),x,y);
end;


procedure C3_to_P(x,y,z:real; var r,l,f:real);
  {cartesian 3-d system to polar one (radians)}
begin
if x=0 then if y>=0 then l:=pi/2 else l:=-pi/2
else
 begin
  l:=arctan(y/x);
  if  x<0  then l:=l+pi
  else if (y<0) then l:=l+2*pi
 end;
r:=sqr(x)+sqr(y);   {here just as an auxiliary real variable}
if r=0 then
 if z>=0 then
  f:=pi/2
 else f:=-pi/2
else
 f:=arctan(z/sqrt(r));
r:=sqrt(r+sqr(z))   {this is r at last}
end;

procedure C3_to_Pd(x,y,z:real; var r,l,f:real);
  {cartesian 3-d system to polar one (degrees)}
begin
 C3_to_P(x,y,z,r,l,f);
 l:=rad_to_deg(l);
 f:=rad_to_deg(f);
end;

procedure P3_to_C(r,l,f:real; var x,y,z:real);
  {polar 3-d system (degrees) to cartesian one}
begin
 z:=cos(f);
 x:=r*cos(l)*z;
 y:=r*sin(l)*z;
 z:=r*sin(f)
end;

procedure P3d_to_C(r,l,f:real; var x,y,z:real);
  {polar 3-d system (degrees) to cartesian one}
begin
 P3_to_C(r,deg_to_rad(l),deg_to_rad(f),x,y,z);
end;

procedure Vector_P3s(V:V3; var l,f:real);
  {3-vector to polar sky}
var r:real;
begin
 C3_to_P(V[1],V[2],V[3],r,l,f);
  {last line of that procedure is obsolete, but that's just one sqrt}
end;

procedure P3s_Vector(L,F:real; {like RA and Decl}
                 var V:V3 {V[3] to max F, V[1] to 0 L});
var C: real;
begin
 V[3]:=sin(F);
 C:=cos(F);
 V[1]:=C*cos(L);
 V[2]:=C*sin(L)
end;

procedure mat(var M:M3 {matrix}; coord:c3d; omega:real{in radians});
          {transformation matrix
           into new frame, that is rotated
            by Omega around Coord}
var C,S:real;
    Cf,Cs:c3d;
begin
 for i:=1 to 3 do for j:=1 to 3 do M[i,j]:=0;
 C:=cos(omega);
 S:=sin(omega);
 Cf:=1+ (coord mod 3);
 Cs:=1+ ((coord+1) mod 3);
 M[coord,coord]:=1;
 M[Cf,Cf]:=C;
 M[Cf,Cs]:=S;
 M[Cs,Cf]:=-S;
 M[Cs,Cs]:=C;
end;

procedure mat_vect(M:M3; var V:V3 {vector});
          {vector:= matrix * vector}
var VA:array[1..3] of real;
begin
 for i:=1 to 3 do VA[i]:=0;
 for i:=1 to 3 do
  for j:=1 to 3 do
   VA[i]:=VA[i]+M[i,j]*V[j];
 for i:=1 to 3 do V[i]:=VA[i]
end;

procedure matl_mat(ML:M3; var M:M3);
          {matrix:= matrix_left * matrix}
var MA: M3;
begin
 for i:=1 to 3 do for j:= 1 to 3 do
  begin
   MA[i,j]:=0;
   for k:=1 to 3 do MA[i,j]:=MA[i,j]+ML[i,k]*M[k,j]
  end;
 for i:=1 to 3 do for j:= 1 to 3 do M[i,j]:=MA[i,j]
end;

procedure matV(var M:M3 {matrix}; VNP_3,V0L_1:V3);
          {transformation matrix into new base}
var l,b: real; MA:M3;
begin
 Vector_P3s(VNP_3,l,b);
 mat(M,3,l);
 mat(MA,2,pi/2-b);
 matl_mat(MA,M);
 mat_vect(M,V0L_1);
 Vector_P3s(V0L_1,l,b);
 mat(MA,3,l);
 matl_mat(MA,M);
end;

procedure matVinv(var M:M3 {matrix}; VNP_3,V0L_1:V3);
          {transformation matrix from the new base to the standard one}
begin
 matV(M,VNP_3,V0L_1);
 VNP_3:=V_unit_base[3];
 V0L_1:=V_unit_base[1];
 mat_vect(M,VNP_3);
 mat_vect(M,V0L_1);
 matV(M,VNP_3,V0L_1);
end;



function  Gnom3 (V:V3 {vector};
                 R:real {distance to the projection plane (along V[3])};
                 var X,Y: real {gnomonic coordinates} )
               : boolean {true if V[3]>0};
begin
 Gnom3:=false;
 if (abs(V[3]))<0.0003 {0.03} then exit;
 X:=R*V[1]/V[3];
 Y:=R*V[2]/V[3];
 Gnom3:= V[3]>0;
end;

procedure Ah_to_td(Pole_rot:real;var L,F:real);
                      {lefthand systems like horizontal->equatorial
                    (then Pole_rot=-(pi/2-Fi), backwards pi/2-Fi), radians}
begin
 mat(Ma_Ah_to_td,2,Pole_rot);
          {transformation matrix for Rotation by Omega around Coord}
 Ah2td(l,f);
end;

procedure Ah2td(var L,F:real); {like ah_to_td, but uses its Ma_...}
var Ve:V3;
begin
 P3s_Vector(-L,F,Ve); {V[3] to max F, V[1] to 0 L}
 mat_vect(Ma_Ah_to_td,Ve);
          {vector:= matrix * vector}
 Vector_P3s(Ve,L,F);
 L:=2*pi-L;
end;

procedure ra_d2l_b(Pole_rot:real;var L,F:real);
                  {righthand systems equatorial->ecliptical-like, radians,
                   pole_rot = epsi_rad for really eq.->ecl. transformation,
                             -epsi_rad for        ecl.->equatorial transf.}
var Ve:V3; Ma:M3;
begin
 P3s_Vector(L,F,Ve); {V[3] to max F, V[1] to 0 L}
 mat(Ma,1,Pole_rot);
          {transformation matrix for Rotation by Omega around Coord}
 mat_vect(Ma,Ve);
          {vector:= matrix * vector}
 Vector_P3s(Ve,L,F);
end;

function object_phase_angleD(l_helioc,b_helioc,l_geoc,b_geoc:real):real;
  {ecliptical coordinates of the object from the Sun and the Earth
   output and input in degrees, ouput is angle Sun-object-Earth}
var pa,ad:real;
begin
 pa:=deg_to_rad(l_geoc-l_helioc);
 ad:=deg_to_rad(b_geoc);
 ah_to_td(deg_to_rad(90-b_helioc),pa,ad);
 object_phase_angleD:=90-rad_to_deg(ad);
end;

procedure ang_dist_posi_angD2(ra1_d,de1_d,ra0_d,de0_d,ra2_d,de2_d:real;var pa,ad: real);
begin
  {all angles convert to radians:}
 pa:=deg_to_rad(-(ra2_d-ra1_d));
 ad:=deg_to_rad(de2_d);
 ah_to_td(deg_to_rad(90-de1_d),pa,ad);
 ra0_d:=deg_to_rad(-(ra0_d-ra1_d));
 de0_d:=deg_to_rad(de0_d);
 ah_to_td(deg_to_rad(90-de1_d),ra0_d,de0_d);
 pa:=Bas_Int(rad_to_deg(pa-ra0_d),360);
 ad:=90-rad_to_deg(ad);
end;

procedure ang_dist(ra1_d,de1_d,ra2_d,de2_d:real;var ad: real);
var pa:real;
begin
  {all angles convert to radians:}
 pa:=deg_to_rad(-(ra2_d-ra1_d));
 ad:=deg_to_rad(de2_d);
 ah_to_td(deg_to_rad(90-de1_d),pa,ad);
 ad:=90-rad_to_deg(ad);
end;

procedure vect_prod_dirD(l1,b1,l2,b2:real; var lp,bp:real);
              {just direction of vector product, all angles in degrees}
var V1,v2,vp:V3;
begin
P3s_vector(deg_to_rad(l1),deg_to_rad(b1),v1);
P3s_vector(deg_to_rad(l2),deg_to_rad(b2),v2);
vect_prod(v1,v2,vp);
vector_p3s(vp,lp,bp);
lp:=rad_to_deg(lp);
bp:=rad_to_deg(bp);
end;

procedure vect_prod_dir(l1,b1,l2,b2:real; var lp,bp:real);
              {just direction of vector product}
var V1,v2,vp:V3;
begin
P3s_vector(l1,b1,v1);
P3s_vector(l2,b2,v2);
vect_prod(v1,v2,vp);
vector_p3s(vp,lp,bp);
lp:=lp;
bp:=bp;
end;

end.