#!/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_HII.".$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.HII.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.HII.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";
}
#
# HII regions
#


$call="correct_vel_map_cube_rot.pl ".$FILE." ve.".$NAME.".vel_map.fits.gz ".$NAME.".V500.zcube.fits.gz ".$vel;
system($call);

$call="correct_vel_map_cube_rot.pl  GAS.".$NAME.".cube.fits.gz ve.".$NAME.".vel_map.fits.gz GAS.".$NAME.".zcube.fits.gz ".$vel;
system($call);
#exit;

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

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

$call="HII_recover.pl map_Ha.".$NAME.".fits.gz 0.4 3.5 0.1 0.1  seg_Ha.".$NAME.".fits.gz mask_Ha.".$NAME.".fits.gz > HII_recover.".$NAME.".log";
mycall($call);
$call="./plot_seg_map.py ".$NAME;
mycall($call);
$call="spec_extract_cube.pl ".$FILE_z." seg_Ha.".$NAME.".fits.gz HII.".$NAME.".rss.fits.gz";
mycall($call);
$call="spec_extract_cube_error.pl ".$FILE." seg_Ha.".$NAME.".fits.gz e_HII.".$NAME.".rss.fits.gz";
mycall($call);

#$d_red=0;

$min_red=$red-20/$vel_light;
$max_red=$red+20/$vel_light;
$d_red=5/$vel_light;

$call="rm -f elines_auto_ssp.HII.".$NAME.".rss.out";
mycall($call);
$call="auto_ssp_elines_several_Av_log_rss_new_error.pl  HII.".$NAME.".rss.fits.gz ".$temp_2." auto_ssp.HII.".$NAME.".rss.out ".$DIR_CONF."/mask_elines.txt ".$DIR_CONF."auto_ssp_V500_several.HII.config ".$plot_rss." -2 5 3850 6800 ".$DIR_CONF."/emission_lines.txt  ".$red." ".$d_red." ".$min_red." ".$max_red." 4.5 0.5 1.2 9.0 0.5 0.2 0.0 1.6";
mycall($call);

$call="gzip -f *.fits";
mycall($call);

$call="plot_output_auto_ssp_elines_several_Av_log_rss_all.pl output.auto_ssp.HII.".$NAME.".rss.out.rss.fits.gz plot_fit.HII.".$NAME.".ps/CPS ";
mycall($call);

$call="get_index_output_auto_ssp_elines_several_Av_log_rss_outFIT3D.pl output.auto_ssp.HII.".$NAME.".rss.out.rss.fits.gz 5 auto_ssp.HII.".$NAME.".rss.out /null > indices.HII.".$NAME.".rss.out";
mycall($call);


$call="cp HII.".$NAME.".rss.pt.txt.gz  HII.".$NAME.".rss.pt.txt";
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";
}