#!/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_2.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_2.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.".$NAME.".log";
open(LOG,">$logfile");

print_time();

$call="get_slice.pl ".$FILE." img ".$DIR_CONF."/slice_V.conf"; 
mycall($call);   
#$call="gzip -d img_V_4500_5500.fits.gz";
#mycall($call);                                                               
$img=rfits("img_V_4500_5500.fits");                                             
($nx,$ny)=$img->dims;                                                           
$val_max=-1e12;                                                                 
$XC=0;                                                                         
$YC=0;                        
$SUM=0;                                                 
for ($ii=30;$ii<45;$ii++) {                                                     
    for ($jj=30;$jj<45;$jj++) {                                                 
	$val=$img->at($ii,$jj);                                                 
	$XC=$XC+$ii*($val)**5;                                                            
	$YC=$YC+$jj*($val)**5;                                                           
	$SUM=$SUM+($val)**5;
    }                                                                           
}   
$XC=$XC/$SUM;
$YC=$YC/$SUM;
$call="cp img_V_4500_5500.fits ".$NAME.".V.fits";
mycall($call);
$call="clean_nan.pl  ".$NAME.".V.fits -1";
mycall($call);
$call="./radial_sum_cube_e.pl ".$FILE." 2.5 ".$XC." ".$YC." rad.".$NAME.".rss.fits 2"; 
mycall($call);                                                                     
$call="./img2spec_e.pl rad.".$NAME.".rss.fits 0 ".$NAME.".spec_5.txt";
mycall($call);                                                                     
$call="./img2spec_e.pl rad.".$NAME.".rss.fits 5 ".$NAME.".spec_30.txt";
mycall($call);                                             

$call="ncftp < ftp.mice";
mycall($call);

#exit;

$call="read_img_header.pl ".$FILE." MED_VEL > SYS_VEL_NOW";
#/disk-b/sanchez/ppak/legacy//DATA/V500/reduced_v1.3c//NGC5947.V500.rscube.fits.gz   5937.49999999998
#$call="dump_header.pl ".$FILE." | grep PIPE_SYSVEL_MEDIAN > SYS_VEL_NOW";
mycall($call);
open(VEL,"<SYS_VEL_NOW");
$dat=<VEL>; chop($dat);
close(VEL);
#($key,$vel)=split(/\|/,$dat);
($key,$vel)=split(" ",$dat);
$vel =~ s/ //g;

$vel=6653.50716406800;

if ($vel>0) {
    $red=$vel/$vel_light;
    $min_red=$red-700/$vel_light;
    $max_red=$red+700/$vel_light;
    $d_red=30/$vel_light;
} else {
    $red=0.02;
    $min_red=0.0045;
    $max_red=0.032;
    $d_red=30/$vel_light;
}

#$red=0.022045;
#$vel=6609.5;
$red=$vel/$vel_light;
$min_red=$red-300/$vel_light;
$max_red=$red+300/$vel_light;
$d_red=30/$vel_light;

$call="auto_ssp_elines_several_Av_log_new.pl  ".$NAME.".spec_5.txt ".$template." auto_ssp.".$NAME.".cen.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.01 0.3 0.0 2.1";
mycall($call);

if ($vel==0) {
    $outfile="auto_ssp.".$NAME.".cen.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";
}



$call="auto_ssp_elines_several_Av_log_new.pl  ".$NAME.".spec_30.txt ".$template." auto_ssp.".$NAME.".int.out ".$DIR_CONF."/mask_elines.txt ".$DIR_CONF."auto_ssp_V500_several.config ".$plot." -3 150 3800 6800 ".$DIR_CONF."/emission_lines.txt  ".$red." ".$d_red." ".$min_red." ".$max_red." 4.5 0.5 1.2 9.0 0.01 0.3 0.0 2.1";
mycall($call);

# REMOVE THIS!
#$WREF=6562.68;
#$WMIN=$WREF*(1+$min_red)-50;
#$WMAX=$WREF*(1+$max_red)+50;
#$call="vel_eline_cube.pl GAS.".$NAME.".cube.fits 1 0.95 ve.".$NAME." ".$WREF." /null ".$WMIN." ".$WMAX;
#mycall($call);

#exit;
#exit;
#exit;

$V_file=$NAME.".V.fits";
$pdl_V=rfits($V_file);
($nx,$ny)=$pdl_V->dims;

#$call="cont_seg_all.pl ".$NAME.".V.fits 0.04 7 0.7 0.02 cont_seg.".$NAME.".fits DMASK.".$NAME.".fits";
#$call="cont_seg_all.pl ".$NAME.".V.fits 0.02 7 0.85 0.005 cont_seg.".$NAME.".fits DMASK.".$NAME.".fits";
#$call="cont_seg_all.pl ".$NAME.".V.fits 0.01 7 0.85 0.001 cont_seg.".$NAME.".fits DMASK.".$NAME.".fits";
#$call="cont_seg_all.pl ".$NAME.".V.fits 0.015 5 0.85 0.00 cont_seg.".$NAME.".fits DMASK.".$NAME.".fits";
#$call="cont_seg_all.pl ".$NAME.".V.fits 0.015 6 0.8 0.00 cont_seg.".$NAME.".fits DMASK.".$NAME.".fits";
$call="cont_seg_all.pl ".$NAME.".V.fits 0.02 7 0.8 0.005 cont_seg.".$NAME.".fits DMASK.".$NAME.".fits";
mycall($call);
#0.02 5 0.85 0.001
$call="spec_extract_cube_mean.pl ".$FILE." cont_seg.".$NAME.".fits CS.".$NAME.".RSS.fits";
mycall($call);

$call="spec_extract_cube_error.pl ".$FILE." cont_seg.".$NAME.".fits e_CS.".$NAME.".RSS.fits";
mycall($call);

$call="rss_seg2cube.pl CS.".$NAME.".RSS.fits cont_seg.".$NAME.".fits SEG.cube.fits";
mycall($call);

$call="get_slice.pl SEG.cube.fits SEG_img ".$DIR_CONF."/slice_V.conf"; 
mycall($call);                                                                  
$call="cp SEG_img_V_4500_5500.fits SEG.V.fits";
mycall($call);                                             
$call="clean_nan.pl  SEG.V.fits -1";
mycall($call);
$call="imarith.pl ".$NAME.".V.fits '/' SEG.V.fits scale.seg.fits";
mycall($call);

print " Fitting RSS file\n";

$call="rm -f elines_auto_ssp.CS.".$NAME.".rss.out";
mycall($call);
$call="auto_ssp_elines_several_Av_log_rss_new_error.pl  CS.".$NAME.".RSS.fits ".$template." auto_ssp.CS.".$NAME.".rss.out ".$DIR_CONF."/mask_elines.txt ".$DIR_CONF."auto_ssp_V500_several.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.3 0.3 0.0 2.1";
mycall($call);
$call="get_index_output_auto_ssp_elines_several_Av_log_rss_outFIT3D.pl output.auto_ssp.CS.".$NAME.".rss.out.rss.fits 5 auto_ssp.CS.".$NAME.".rss.out /null > indices.CS.".$NAME.".rss.out &";
mycall($call);


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


$call="./index_seg_cube.pl indices.CS.".$NAME.".rss.out cont_seg.".$NAME.".fits indices.CS.".$NAME.".cube.fits &";
mycall($call);

#print "DONE\n";

$call="map_auto_ssp.seg.pl elines_auto_ssp.CS.".$NAME.".rss.out cont_seg.".$NAME.".fits map.CS.".$NAME." &";
mycall($call);

$call="FIT3D_output_rss_seg2cube.pl output.auto_ssp.CS.".$NAME.".rss.out.rss.fits 1 cont_seg.".$NAME.".fits SSP_mod.cube.fits";
mycall($call);

$call="imarith.pl SSP_mod.cube.fits '*' scale.seg.fits SSP_mod.".$NAME.".cube.fits";
mycall($call);

$call="imarith.pl ".$FILE." '-' SSP_mod.".$NAME.".cube.fits TMP.".$NAME.".cube.fits";
mycall($call);
$call="rm -r test.cube.fits";
mycall($call);
$call="smooth_spec_clip_cube.pl TMP.".$NAME.".cube.fits test.cube.fits 75 1.5 10 1860";
mycall($call);
$call="imarith.pl TMP.".$NAME.".cube.fits - test.cube.fits GAS.".$NAME.".cube.fits";
mycall($call);





print " GAS.$NAME.cube.fits created\n";

$WREF=6562.68;
$WMIN=$WREF*(1+$min_red)-50;
$WMAX=$WREF*(1+$max_red)+50;
$call="vel_eline_cube.pl GAS.".$NAME.".cube.fits 1 0.95 ve.".$NAME." ".$WREF." /null ".$WMIN." ".$WMAX;
mycall($call);

#
# Analyze cube
#
$FILE_CONF_NOW=$DIR_CONF."/auto_ssp_V500_several_SII.config";
open(CONF,"<$FILE_CONF_NOW");
$conf=<CONF>;
$conf=<CONF>;
$conf=<CONF>;
$nr=<CONF>;
chop($nr);
print "NR=$nr\n";
for ($j=0;$j<$nr;$j++) {
    $conf=<CONF>;
    chop($conf);
    @d_conf=split(" ",$conf);
    $w1=$d_conf[0];
    $w2=$d_conf[1];
    $in_conf=$d_conf[3];
    open(IN,"<$in_conf");
    open(TMP,">tmp.conf");
    $ni=0;
    while($in=<IN>) {
	chop($in);
	if ($ni!=5) {
	    print TMP "$in\n";
	} else {
	    $vel1=$vel-300;
	    $vel2=$vel+300;
	    print TMP "$vel     1       $vel1     $vel2   -1\n";
	}
	$ni++;
    }
    close(TMP);
    close(IN);
# ve.NGC5947.mask_map.fits
#    $call="kin_back_cube.pl  GAS.".$NAME.".cube.fits tmp.conf none 0 ".$w1." ".$w2." KIN.GAS.".$w1."_".$w2.".".$NAME.".out /null none 1 > junk.fit.log ";
    if ($j==0) {
	$call="./kin_back_cube_guided_error.pl  GAS.".$NAME.".cube.fits ".$FILE." tmp.conf none 0 ".$w1." ".$w2." KIN.GAS.".$w1."_".$w2.".".$NAME.".out /null  ve.".$NAME.".vel_map.fits ve.".$NAME.".mask_map.fits 1 > junk.fit.log ";
    } else {
#	$call="./kin_back_cube_fixed_error.pl  GAS.".$NAME.".cube.fits ".$FILE." tmp.conf none 0 ".$w1." ".$w2." KIN.GAS.".$w1."_".$w2.".".$NAME.".out /null  map.6530_6810.".$NAME."_vel_00.fits disp.fits 1 > junk.fit.log ";
	$call="./kin_back_cube_fixed_error.pl  GAS.".$NAME.".cube.fits ".$FILE." tmp.conf none 0 ".$w1." ".$w2." KIN.GAS.".$w1."_".$w2.".".$NAME.".out /null  vel.fits disp.fits 1 > junk.fit.log ";
    }
    mycall($call);
#map.6530_6950.NGC5947_vel_00.fits

    $call="mapgrid_back_rot.pl ".$nx." ".$ny." KIN.GAS.".$w1."_".$w2.".".$NAME.".out map.".$w1."_".$w2.".".$NAME." 1 > junk.fit.log";
    mycall($call);
    if ($j==0) {
	$call="rm -f vel.fits";
	mycall($call);
	$call="rm -f disp.fits";
	mycall($call);
	$call="med2df.pl map.6530_6810.".$NAME."_vel_00.fits vel.fits 3 3";
	mycall($call);
	$call="med2df.pl map.6530_6810.".$NAME."_disp_00.fits disp.fits 4 4";
	mycall($call);
    }

    print "$w1-$w2 region analyzed\n";
}
close(CONF);


#
# HII regions
#

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

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


$call="spec_extract_diffuse.pl  ".$FILE."  map_Ha.".$NAME.".fits diffuse.fits diff.".$NAME.".rss.fits 0.01";
mycall($call);
$call="cp seg.diffuse.fits seg_diff.".$NAME.".fits";
mycall($call);

$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 ".$temp_2." auto_ssp.HII.".$NAME.".rss.out ".$DIR_CONF."/mask_elines.txt ".$DIR_CONF."auto_ssp_V500_several.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.3 0.3 0.0 2.1";
mycall($call);
$call="plot_output_auto_ssp_elines_several_Av_log_rss_all.pl output.auto_ssp.HII.".$NAME.".rss.out.rss.fits 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 5 auto_ssp.HII.".$NAME.".rss.out /null > indices.HII.".$NAME.".rss.out";
mycall($call);

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

$call="ncftp < ftp.mice";
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";
}