#!/usr/bin/perl
#
#
# This program find peaks in a 2D fiber based spectral image
#
#

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;
$ENV{PGPLOT_FOREGROUND} = "black";
$ENV{PGPLOT_BACKGROUND} = "white";


require "/home/sanchez/sda1/perl/MY/my.pl";

$galfit="/home/sanchez/sda1/galfit/galfit";
$cl="/home/sanchez/sda1/iraf/iraf/unix/hlib/cl.csh";

if ($#ARGV<2) {
    print "USE: plot_output_auto_ssp_elines_several_Av_log_rss.pl INPUT_FILE.FITS F_CUT wavelength DEVICE [MIN MAX] [WMIN WMAX]\n [REF_LINES CUT]";
    exit;
}

$input=$ARGV[0];
$f_cut=$ARGV[1];
$w_cut=$ARGV[2];
$dev=$ARGV[3];
$def=0;
if ($#ARGV==5) {
    $min=$ARGV[4];
    $max=$ARGV[5];
    $def=1;
}





$y_min=1e12;
$y_max=-1e12;
$n=0;
$pdl=rfits("$input");
($nx,$ny,$nz)=$pdl->dims;
#print "($nx,$ny,$nz)\n";
$crval=$pdl->hdr->{"CRVAL3"};
$cdelt=$pdl->hdr->{"CDELT3"};
$crpix=$pdl->hdr->{"CRPIX3"};





$nx0=int(($w_cut-50-$crval)/$cdelt);
$nx1=int(($w_cut+50-$crval)/$cdelt);
$pdl_cut=$pdl->slice("$nx0:$nx1,:,(0)");


$flux = medover($pdl_cut);
@b=$flux->dims;

$NY=0;
while((($flux->at($NY)>$f_cut))&&($NY<($b[0]-1))) {
    $val=$flux->at($NY);
    $NY++;
}

$pdl_org=$pdl->slice(":,0:$NY,(0)");
$pdl_mod=$pdl->slice(":,0:$NY,(1)");

$pdl_ratio=$pdl_org/$pdl_mod;
$pdl_rot = $pdl_ratio->xchg(0,1);
($mean,$prms,$median,$min,$max,$adev,$rms) = statsover($pdl_rot);

@c=$mean->dims;

#print "B=@b, NY=$NY, C=@c\n";


#print "$crval $cdelt $crpix\n";
if ($cdelt==0) {
    $cdelt=1;
}
open(OUT,">rat_mod.out");
for ($i=0;$i<$nx;$i++) {
    $wave[$i]=$crval+$cdelt*($i+1-$crpix);
    $f_mean[$i]=$mean->at($i);
    $f_median[$i]=$median->at($i);
    $f_min[$i]=$min->at($i);
    $f_max[$i]=$max->at($i);
    $f_rms[$i]=$rms->at($i);
    $f_0[$i]=$f_mean[$i]-$f_rms[$i];
    $f_1[$i]=$f_mean[$i]+$f_rms[$i];
    if ($f_mean[$i]>$y_max) {
	    $y_max=$f_mean[$i];
    }
    if ($f_mean[$i]<$y_min) {
	$y_min=$f_mean[$i];
    }
    print OUT "$i $wave[$i] $f_mean[$i]\n";
}
close(OUT);

$wmin=$wave[0];
$wmax=$wave[$nx-1];
#print "$wmin $wmax\n";
if ($#ARGV==5) {
    $min=$ARGV[4];
    $max=$ARGV[5];
    $def=1;
}

if ($#ARGV==7) {
    $min=$ARGV[4];
    $max=$ARGV[5];
    $wmin=$ARGV[6];
    $wmax=$ARGV[7];
    $def=1;
}
$ref_def=0;
if ($#ARGV==8) {
    $min=$ARGV[4];
    $max=$ARGV[4];
    $wmin=$ARGV[6];
    $wmax=$ARGV[7];
    $LABEL=$ARGV[8];
    $def=1;
}

if ($def==1) {
    $y_min=$min;    
    $y_max=$max;
}

#print "$y_min $y_max\n";

pgbegin(0,$dev,1,1);
pgsfs(1.2);
pgscf(2);             # Set character font
pgslw(2);             # Set line width
pgsch(1.2);           # Set character height
pgenv($wmin,$wmax,$y_min,$y_max,0,0);
pglabel("Wavelength (\\A)","F\\dorg\\u/F\\dmod\\u","$input");
pgsci(1);
pgsls(1);
pgline($nx,\@wave,\@f_median);
pgsls(2);
pgline($nx,\@wave,\@f_mean);
pgsci(8);
pgsls(2);
pgline($nx,\@wave,\@f_0);
pgsls(2);
pgline($nx,\@wave,\@f_1);


my @sec;
$nn=0;
for ($i=0;$i<$nx;$i++) {
    if (($wave[$i]>$wmin)&&($wave[$i]<$wmax)) {
	$sec[$nn]=$f_mean[$i];
	$nn++;
    }
}


$median=median(@sec);
$mean=mean(@sec);
$sig=sigma(@sec);

$flux=$mean*$n*($wave[1]-$wave[0]);

print "$median +- $sig (F=$flux)\n";
pgsls(1);
pgptxt($wmin+0.05*($wmax-$wmin),$ymax-0.05*($ymax-$ymin),0,0,"$mean+-$sig\n");
pgclose;
pgend;




exit;


sub PLOT {
    for ($i=0;$i<$nx;$i++) {
	$wave_now[$i]=$wave[$i][$k];
	$flux_now[$i]=$flux[$i][$k];
    }
    $color=$k+1;
    pgsci($color);
    if ($color==1) {
	pgslw(3);
    } else {
	pgslw(2);
    }
    pgline($nx,\@wave_now,\@flux_now);
#    $color++;
#    if ($color>15) {
#	$color=1;
#    }
}