#!/usr/bin/perl

# konstante observatorija
# OBS=('46.04361','14.5272')
# AGO
$Longitude=0.9684889;
$Latitude=46.043056;




# za pretvorbo st v radijane
$rad=0.0174532925199433;

$home = $ENV{"HOME"};     #  "Home" imenik
#$home    ="/home/bojan/tmp/srctmp/pika/preglej/acit";

$bin=2;
$rep=1;
#--------------------------------------------------------#
# seznam filtrov za pripravo request datoteke            #
#                                                        #
$filter[0] = 'R';                                        #
$filter[1] = 'V';                                        #
$filter[2] = 'B';                                        #
$filter[3] = 'I';                                        #
#                                                        #
#--------------------------------------------------------#

#--------------------------------------------------------#
# filter in ekspozicija za match sliko
#
$match_filter = 'V';
$match_expo   = 10;
#--------------------------------------------------------#

#--------------------------------------------------------#
# ekspozicijski casi za posamezni filter                 #
#                                                        #
$exp   [0] = 15;                                         #
$exp   [1] = 20;                                         #
$exp   [2] = 60;                                         #
$exp   [3] = 20;                                         #
#                                                        #
#--------------------------------------------------------#




$EW = 1.0;
if ($ARGV[0] eq "-w"){
	$EW = -1.0;
	shift;
}	
if ($ARGV[0] eq "-e"){
	shift;
}	
$rqsfile = $ARGV[0];


$mjd = MJD();
if($rqsfile eq ""){
	$rqsfile=sprintf "landolt-%9.3f",$mjd;
}	
print "request datoteka:$rqsfile\n";


#open (LF,"$home/bin/landolt-fields.dat");
open (LF,"$home/bin/landolt-polja-2015.dat");

print "polje                    ra             dec              X            H         visina\n";
      #SA100                  08:53:36       -00:43:26       8.8567         5.37         5.98

# format bi naj bil
# SA92                   00:56:11       +00:42:52   V    2   30   1 

$i = 0;
while(<LF>){
    if(/^#/) {next}	
	# prebere datoteko in preracuna zracne mase zdaj
	# format datoteke:
	# SA92 00:56:11 +00:42:52 2000.0
	# razsirjeni format, dodamo koordinate centra in cas ekspozicije ( naredi V posnetek)
	# SA93  01:54:36 +00:39:24 2000.0 8.83  01:50:36 +00:50:24   1
	# datoteka 2009
	# SA92 		00:55:15.5 +00:47:26 70x70 10.60   15.92  +0.32   +1.64 

#	($ime[$i],$ra_s[$i],$dec_s[$i],$epo[$i],$Vmag[$i],$ra_center_s[$i],$dec_center_s[$i],$exp_center[$i])=split;
	($ime[$i],$ra_s[$i],$dec_s[$i],$polje,$mag1,$mag2,$bv1,$bv2)=split;

	$ra[$i] = strtora($ra_s[$i])/3600;
	$dec[$i]= strtodec($dec_s[$i])/3600;
	$h[$i]  = Alt($mjd,$ra[$i],$dec[$i]);
	if($h[$i]>3){
		# pozor formula za izracun velja za X < 6.7
		$X[$i]  = airmass($mjd,$ra[$i],$dec[$i]);
	}
	$H[$i]  = $Lst - $ra[$i]/15;  # /

	if($H[$i]>12){$H[$i]-=24}
	if($H[$i]<-12){$H[$i]+=24}
	if($X[$i]>0){	
		printf "%-15s %15s %15s %12.4f %12.2f %12.2f\n",$ime[$i],$ra_s[$i],$dec_s[$i],$X[$i],$H[$i],$h[$i];
	}
	$i++;
}
close (LF);

#$n = select_fields_bin();
#if ($n <3){
$n = select_fields_minmax();
#}	

#--------------------------------------------------#
# zapise datoreko request                          #
#                                                  #
#--------------------------------------------------#
open(RS,">snemalna-landolt.txt");

#
# izbrane zvezde so oznacene 1,2,3
for($j=$n;$j>0;$j--){
	
	for($i=0;$i<=$#ime;$i++){	

		if($N[$i] eq $j){
        	#
        	# najprej naredi match in popravi pozicijo teleskopa
        	#
          #  printf RS "MATCH      %12s %12s %6.0f %5.1f %1s\n",$ra_s[$i],$dec_s[$i],$epo[$i],$match_expo,$match_filter;        	
        	# center na svetlo zvezdo, opusceno
        	#if(length($ra_center_s[$i])>4 && length($dec_center_s[$i])>4){
        		# ce so vnesene koordinate centra
        	#	printf RS "CENTER     %12s %12s %6.0f %5.2f V\n",$ra_center_s[$i],$dec_center_s[$i],$epo[$i],$exp_center[$i];
		    #}
			for($f=0;$f<=$#filter;$f++) {
				# za vse filtre
	     		#printf RS "%-10s %12s %12s %6.0f %5.1f %1s\n",$ime[$i]."_".$filter[$f],$ra_s[$i],$dec_s[$i],$epo[$i],$exp[$f],$filter[$f];
					printf RS "%-10s %12s %12s %1s %1d %3d %3d \n",$ime[$i],$ra_s[$i],$dec_s[$i], $filter[$f],$bin,$exp[$f],$rep;
			} 
			last;
      }		
	}
}
close(RS);

exit;






sub select_fields_minmax {
#-----------------------------------------------------#
# izberi 3 polja po metodi minmax                     #
# izberi polje najvisje,najnizje v sredini            #
# boljse prvi in zadnja dva						      #
#-----------------------------------------------------#

print "\nizbira:\n";

$x=1;
$n=1;
# uredimo po narascajocem X
while ($x < 4.4){
	for($i=0;$i<=$#ime;$i++){
		if($ix[$i] < 1){	
			if($EW * $H[$i]<0){
				if( $X[$i] > $x &&  $X[$i] <  $x+0.001){
					$ix[$i] = $n;					
					print " $n     $X[$i] \n";
					$n++;
					last;
				}	
			}
		}
	}
	$x += 0.001;
}	
for($i=0;$i<=$#ime;$i++){
	if($ix[$i] == $n-1){
		# zadnji je najvecji
		$N[$i] = 3;
	}
	if($ix[$i] == $n-2){
		# predzadnji, srednji pa bi bil int($n/2)
		$N[$i] = 2;
	}
	if($ix[$i] == 1){
		# prvi
		$N[$i] = 1;
	}	
	if($N[$i]>0){
		printf "%-15s %15s %15s  %12.4f %12.2f %12.2f \n",$ime[$i],$ra_s[$i],$dec_s[$i],$X[$i],$H[$i],$h[$i];
	}
}
return 3;
}


sub select_fields_bin {
#-----------------------------------------------------#
# izberi 3 polja po metodi 3 bining                   #
# idealno 1.4 2.8 4.2                                 #
# vrne stevilo izbranih polj						  #
# visina X=4 -> 14 st                                 #
#        X=5 -> 12 st								  #
#-----------------------------------------------------#
$meja[0]  = 1.4;    # 1-1.8        0.8
$meja[1]  = 2.5;    # 1.8 - 3.2    1.4
$meja[2]  = 3.5;    # 2.7 - 4.3    1.6

$eps [0] = $eps[1] = $eps[2]  = 0;
$n       = 1;     # stevilo nabranih zvezd

print "\nizbira:\n";
while($n<4 && $eps[0]<0.4){
	
	for($i=0;$i<=$#ime;$i++){	
		
		if($EW * $H[$i]<0){
			# vzhodna ali zahodna staran
			for($j=0;$j<3;$j++){				
				if($N[$i] < 1){
					# ta zvezda se ni bila izbrana
					if($X[$i]>$meja[$j]-$eps[$j] && $X[$i]<$meja[$j]+$eps[$j]){
						printf "%-15s %15s %15s %12.4f %12.2f %12.2f \n",$ime[$i],$ra_s[$i],$dec_s[$i],$X[$i],$H[$i],$h[$i];	
						$N[$i] = $n;			
						$meja[$j] = -1e9;
						$n++;
						last;
					}
				}
    		}
	    }
	}
    $eps [0]  += 0.008;
	$eps [1]  += 0.014;
	$eps [2]  += 0.014;
}
print "izbrano stevilo polj:$n\n";
return $n-1;
}










sub asin { atan2($_[0], sqrt(1 - $_[0] * $_[0])) }
sub acos { atan2( sqrt(1 - $_[0] * $_[0]), $_[0] ) }

sub MJD
#---------------------------------------*/
# MJD                                   */
# ura je decimalni del dneva            */
#---------------------------------------*/
{
	($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = gmtime(time);
	$MJD=(-678987.0+367.0*(1900+$year)+int(int(275*($mon+1)/9)+$mday-int(7*(($year+1900)+int((($mon+1)+9)/12))/4)));
	$MJD+=($sec/3600+$min/60+$hour)/24.0;
	return $MJD;
}

sub LST
#---------------------------------------*/
# krajevni zvezdni cas   				*/
#    							        */
#---------------------------------------*/
{
	($MJD)=@_;
	$UT=( $MJD - int($MJD))*24.0;
	$T =( $MJD - 51544.5 )/36525.0;
	$Gmst = 6.697374558+1.0027379093*$UT+(8640184.812866+(0.093104-0.0000062*$T)*$T)*$T/3600.0;
	#Local mean sideral time
	$Lst= ($Gmst+$Longitude);
    $Lst= $Lst-24.0*int($Lst/24.0);
    if($Lst<0.0){
      $Lst+=24.0;
    }
    return $Lst;
}

sub Alt
#---------------------------------------*/
# Visina nad obzorjem                   */
# rek, dec v stopinjah                  */
#---------------------------------------*/
{
	($MJD,$REK,$DEC)=@_;
   $Lst = LST($MJD);
	$slat = sin($Latitude*$rad);
	$clat = cos($Latitude*$rad);
	$z= ($slat*sin($DEC*$rad)+$clat*cos($DEC*$rad)*cos((15*$Lst-$REK)*$rad) );
	$alt = asin($z)/$rad;
	return $alt;
}


sub airmass
#---------------------------------------*/
# Zracna masa                           */
# rek, dec v stopinjah                  */
#---------------------------------------*/
{
	($MJD,$REK,$DEC)=@_;
   $Lst = LST($MJD);
	$slat = sin($Latitude*$rad);
	$clat = cos($Latitude*$rad);
	$z= ($slat*sin($DEC*$rad)+$clat*cos($DEC*$rad)*cos((15*$Lst-$REK)*$rad) );
	$alt = asin($z)/$rad;
	# z je sin(z)
	$X=0.0;
	if ($z>0.02) {
	   $X = (1.0/$z);
		$X = ($X-0.0018167*($X-1.)-0.002875*($X-1.)*($X-1.)-0.0008083*($X-1.)*($X-1.)*($X-1.));
	}
	return $X;
}


sub strtora

{
    my ($ras) = @_;
    my ($h, $m, $s);
    ($h,$m,$s) = split(":",$ras);
    return 15 * ($h * 3600.0 + $m * 60.0 + $s);
}

#  Transform a string to declination in arseconds
sub strtodec

{
    my ($decs) = @_;
    my ($d, $m, $s)=split(":",$decs);
    if(index($decs,"-") > -1){
      return -1.0*(abs($d) * 3600.0 + abs($m) * 60.0 + abs($s));
    }
    else {
       return (abs($d) * 3600.0 + abs($m) * 60.0 + abs($s));
    }

}
