#!/usr/bin/perl
use Statistics::OLS;
use Math::FFT;
use Math::Stat;
use Math::Spline qw(spline linsearch binsearch);
use Math::Derivative qw(Derivative2);
use Math::Approx;
use Astro::FITS::CFITSIO qw( :longnames :constants );
use PDL;
use PDL::Fit::Polynomial; 
use PDL::Filter::Linear;
use PGPLOT;  # Load PGPLOT module
use PDL::Fit::Gaussian;
use PDL::Slatec;
use PDL::Image2D;
#use PDL::Matrix;
$vel_light=299792.458;

$plot=2;
$plot_rss=0;

#$DIR_WORK="/disk-a/sanchez/ppak/legacy/DATA/SSP/SEG/";
#$DIR_CONF="/disk-b/sanchez/ppak/legacy/";
$DIR_CONF="../legacy/";
#/disk-b/sanchez/ppak/legacy/
#$DIR_DATA="/disk-b/sanchez/ppak/legacy/DATA/V500/reduced_v1.3c/";
$DIR_DATA="/disk-a/sanchez/ppak/legacy/DATA/V500/reduced_v1.3c/";
#chdir($DIR_WORK);
$max_size=10;
$frac=0.9;
$template="/disk-b/sanchez/ppak/legacy/miles_12.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_12.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_6.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_12.fits";
$temp_2="/disk-b/sanchez/ppak/legacy/miles_12.fits";
#
#$template="/disk-b/sanchez/ppak/legacy/miles_6.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_20.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_4.fits";

if ($#ARGV<0) {
    print "ana_single.pl NAME\n";
    exit;
}
$NAME=$ARGV[0];
$FILE=$DIR_DATA."/".$NAME.".V500.rscube.fits.gz";

$logfile="ana_single_map_Ha.".$NAME.".log";
open(LOG,">$logfile");

print_time();

$call="read_img_header.pl ".$FILE." MED_VEL > SYS_VEL_NOW";
mycall($call);
#print "$call\n"; exit;
open(VEL,"<SYS_VEL_NOW");
$dat=<VEL>; chop($dat);
close(VEL);
#($key,$vel)=split(/\|/,$dat);
($key,$vel)=split(" ",$dat);
$vel =~ s/ //g;
if ($vel>0) {
    $red=$vel/$vel_light;
    $min_red=$red-300/$vel_light;
    $max_red=$red+300/$vel_light;
    $d_red=30/$vel_light;
} else {
    $red=0.02;
    $min_red=0.0045;
    $max_red=0.032;
    $d_red=30/$vel_light;
}
$call="auto_ssp_elines_several_Av_log_new.pl  ".$NAME.".spec_5.txt ".$template." auto_ssp.".$NAME.".cen.map_Ha.out ".$DIR_CONF."/mask_elines.txt ".$DIR_CONF."auto_ssp_V500_several.config ".$plot." -3 50 3800 6800 ".$DIR_CONF."/emission_lines.txt  ".$red." ".$d_red." ".$min_red." ".$max_red." 4.5 0.5 1.2 9.0 0.5 0.1 0.0 1.6";
mycall($call);

if ($vel==0) {
    $outfile="auto_ssp.".$NAME.".cen.map_Ha.out";
    open(RED,"<$outfile");
    $RED=<RED>; $RED=<RED>; chop($RED);
    close(RED);
    @RED_DATA=split(" ",$RED);
    $red=$RED_DATA[4];
    $min_red=$red-300/$vel_light;
    $max_red=$red+300/$vel_light;
    $d_red=30/$vel_light;
    $vel=$red*$vel_light;
    print "RED=$red VEL=$vel\n";
}
#
# map_Ha regions
#

$FILE_z=$NAME.".V500.zcube.fits.gz";

#$call="get_Ha.pl ".$FILE_z." ".$red." map_Ha.".$NAME.".fits.gz";
$call="get_Ha.pl ".$FILE." ".$red." map_Ha.".$NAME.".fits.gz";
mycall($call);

print_time();
close(LOG);
exit;


sub print_time {
    my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
    print "# TIME $sec $min $hour $mday $mon $year $wday $yday $isdst\n";
    print LOG "# TIME $sec $min $hour $mday $mon $year $wday $yday $isdst\n";
}

sub mycall {                                                                                                     
    my $call=$_[0];
    print "$call | PROCESSING\n";
    system($call);
    print LOG "$call\n";
    print "DONE\n";
}