package aux;
# Copyright (C) 2008  Bernhard W. Adams

# 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 3 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, see <http://www.gnu.org/licenses/>.

use strict;
use POSIX;

#################################
#
#
#################################

sub date_to_epoch {
 my ($year,$month,$day)=@_;
 my ($yr_d,$epoch);

 $yr_d=365.256363051;
 if($month>11) {$day+=30;}
 if($month>10) {$day+=31;}
 if($month>9) {$day+=30;}
 if($month>8) {$day+=31;}
 if($month>7) {$day+=31;}
 if($month>6) {$day+=30;}
 if($month>5) {$day+=31;}
 if($month>4) {$day+=30;}
 if($month>3) {$day+=31;}
 if($month>2) {if(leapyr($year)) {$day+=29;} else {$day+=28;}}
 if($month>1) {$day+=31;}
# (J-2000)*$yr_d+2451545.0=JEP
 $epoch=($year-2000)*$yr_d+2451545.0+$day;
# BEP=(B-1900)*365.242198781+2415020.31352
# $epoch=($year-1900)*365.242198781+2415020.31352+$day;
}

#################################
#
#
#################################
sub epoch_to_date {
 my ($epoch)=@_;
 my ($yr_d,$yr,$yep,$dy,$mo,$res);

 $yr_d=365.256363051;
 $yr=2000+floor(($epoch-2451545.0)/$yr_d);
 $yep=($yr-2000)*$yr_d+2451545.0;
 $dy=$epoch-$yep;
# print "dy = ".$dy."\n";
 $mo=0;
 if($mo==0) {if($dy>31) {$dy-=31;} else {$mo="01";}}
 if(leapyr($yr)) {if($mo==0) {if($dy>29) {$dy-=29;} else {$mo="02";}}}
 else            {if($mo==0) {if($dy>28) {$dy-=28;} else {$mo="02";}}} 
 if($mo==0) {if($dy>31) {$dy-=31;} else {$mo="03";}}
 if($mo==0) {if($dy>30) {$dy-=31;} else {$mo="04";}}
 if($mo==0) {if($dy>31) {$dy-=31;} else {$mo="05";}}
 if($mo==0) {if($dy>30) {$dy-=31;} else {$mo="06";}}
 if($mo==0) {if($dy>31) {$dy-=31;} else {$mo="07";}}
 if($mo==0) {if($dy>31) {$dy-=31;} else {$mo="08";}}
 if($mo==0) {if($dy>30) {$dy-=31;} else {$mo="09";}}
 if($mo==0) {if($dy>31) {$dy-=31;} else {$mo="10";}}
 if($mo==0) {if($dy>30) {$dy-=31;} else {$mo="11";}}
 if($mo==0) {if($dy>31) {$dy-=31;} else {$mo="12";}}
 $dy=int($dy*100)/100;
 if($dy<10) {$dy="0".$dy;}
 $res=$yr."/".$mo."/".$dy;
}

##################
#
#
##################
sub leapyr {
 my $yr=$_[0];
 my $res=0;
 
 if($yr%4 == 0) {$res=1;}
 if($yr%100 == 0) {$res=0;}
 if($yr%400 == 0) {$res=1;}
 $res;
}
#################################
#
#################################
sub sec_to_date {
 my ($sec)=@_;
 epoch_to_date($sec/(3600*24));
}
#################################
#
#################################
sub date_to_sec {
 my (@date)=@_;
 my $res;
 $res=24*3600*date_to_epoch(@date);
}

#######################
#
#
#######################
sub sc_prod {
 my ($x,$y,$z,$X,$Y,$Z)=@_;
 my $res;
 $res=$x*$X+$y*$Y+$z*$Z;
 return $res;
}
#######################
#
#
#######################
sub vsum {
 my ($x,$y,$z,$X,$Y,$Z)=@_;
 my @res;
 @res=($x+$X,$y+$Y,$z+$Z);
 return @res;
}
#######################
#
#
#######################
sub vdiff {
 my ($x,$y,$z,$X,$Y,$Z)=@_;
 my @res;
 @res=($x-$X,$y-$Y,$z-$Z);
 return @res;
}

#######################
#
#
#######################
sub vabs {
 my ($x,$y,$z)=@_;
 my $res;
 $res=sqrt($x*$x+$y*$y+$z*$z);
}
#######################
#
#
#######################
sub v_prod_s {
 my ($x,$y,$z,$s)=@_;
 my @res;
 @res=($x*$s,$y*$s,$z*$s);
}
##############
# 
# vector cross product
# 
##############
sub cross {
 my ($ax,$ay,$az,$bx,$by,$bz)=@_;
 my @res;
  
 $res[0]=$ay*$bz-$by*$az;
 $res[1]=$az*$bx-$bz*$ax;
 $res[2]=$ax*$by-$bx*$ay;
 @res;
}
      
##############
#
# vector rotation a about b
#
##############
sub rotate {
 my ($x,$y,$z,$Ax,$Ay,$Az,$angle,$d)=@_;
 my ($aA,@X,@Y);
 my @res;
    
 $aA=$x*$Ax+$y*$Ay+$z*$Az; # component of (x,y,z) along axis (Ax,Ay,Az)
 # part of (x,y,z) in plane perp. to (Ax,Ay,Az), orth. to (x,y,z) 
 $Y[0]=$Ay*$z-$Az*$y;
 $Y[1]=$Az*$x-$Ax*$z;
 $Y[2]=$Ax*$y-$Ay*$x;
 # part of (x,y,z) in plane perp. to (Ax,Ay,Az), "along" to (x,y,z) 
 $X[0]=$Y[1]*$Az-$Y[2]*$Ay; 
 $X[1]=$Y[2]*$Ax-$Y[0]*$Az; 
 $X[2]=$Y[0]*$Ay-$Y[1]*$Ax;
# if($d==1) {print "angle = ".$angle."\n";}
 $res[0]=$aA*$Ax+$X[0]*cos($angle)+$Y[0]*sin($angle);
 $res[1]=$aA*$Ay+$X[1]*cos($angle)+$Y[1]*sin($angle);
 $res[2]=$aA*$Az+$X[2]*cos($angle)+$Y[2]*sin($angle);
# if($d==1) {print "res = (".$res[0].", ".$res[1].", ".$res[2].")\n";} 
 return @res;
}
	 
################################
sub ems {
 my $ms=$_[0];

 die $ms;
}

sub dumpvars {
 my @in=@_;
 my ($m,$n,$pad,@p);
 
 if($in[0] ne "ARRAY") {
   @p=split(/\s+/,$in[0]);
   for($m=0;$m<scalar(@p);++$m) 
     {$pad=3-length($p[$m]); for($n=0;$n<$pad;++$n) {$p[$m]=$p[$m]." ";}}
   for($m=0;$m<scalar(@p);++$m) {
     print $p[$m]." = ".$in[$m+1]."\n";
   }
 }
 else {
   open F,">".$in[1];
   for($m=2;$m<scalar(@in);++$m) {
#     print "el. ".($m-2)." : ".$in[$m]."\n";
     print F ($m-2)."\t".$in[$m]."\n";   
   }
#   print "\n\n";
   close F;
 }
}

#############################
#
#
#############################
sub msg {
 my $ms=$_[0];
 my $len;
 
 $ms=$ms."\n";
 $len=length($ms);
 syswrite(STDOUT,$ms,$len);
}
#############################
#
#
#############################
sub die {
 my ($pr,$tr)=@_;
 
 die "error in ".$pr." called from ".$tr;
}

1; # return value