#!/usr/bin/perl -w
use PDL;
use Vtools;
use PDL::FFT;
use PDL::NiceSlice;
use Linfit_norm;
use PDL::Graphics::PGPLOT::Window;
use PDL::Transform;
$pi = 2*acos(0);


$im = rfits('sfd_zoph.fits.gz');
#1h40g30x30.fits
$im = rfits('1h40g30x30.fits');

$hdr = $im->gethdr;

$N=$$hdr{'NAXIS1'};
print "N: $N \n";
$dx=(1/($N-1));
$x=sequence($N)*$dx;
$xs=$dx*($im->xvals);
$ys=$dx*($im->yvals);
$x_0=$xs($N/2,$N/2)->sclr; $y_0=$ys($N/2,$N/2)->sclr;
$rad=sqrt(($xs-$x_0)**2+($ys-$y_0)**2);

$im *= exp(-$rad**2/(2.*(70*$dx)**2));

&Vtools::view($im);

$hr = $im -> copy;
$hi = $hr*0;

fftnd($hr,$hi);


$hr=&restorefft($hr);
$hi=&restorefft($hi);



$powerspec=sqrt($hr**2+$hi**2);
$f=(sequence($N)-($N/2))/($N);

&Vtools::view($powerspec); 

$Nbins = 50 ; # $N/2;
$f = zeroes($Nbins);
$P = zeroes($Nbins);
foreach $i (0..$Nbins-1) {
    $f($i) .= ($i+1)*($N/(2*$Nbins))*$dx;
    $P($i) .= average($powerspec->where( ($rad<$f(($i)) &  ($rad< ($f(($i)+$dx))))));
}




$constoff = ones(nelem($f)) ; 
$funcs = cat $constoff,log($f);

$sigs=1/($constoff*0.001*max(log($P)));

($yfit,$coefslin,$covlin) = linfit1d_norm log($P),$funcs,{Weights=>$sigs};

wcols '%5.2e$\pm$%5.2e(%4.1f)', $coefslin,sqrt($covlin->diagonal(0,1)),$coefslin/sqrt($covlin->diagonal(0,1));

#$fracdim = sprintf("%5.2f",$coefslin(1) ->sclr);

$index = $coefslin(1) ->sclr;

$xtext=max($f)/2.;
$ytext=max($P)/2.;
$codstr ='$win->text(\' \ga='.sprintf("%6.2f",$index)."\',$xtext,$ytext,{CHARSIZE=>2})";
print "CODE $codstr \n";
&Vtools::spec($f,$P,{code=>$codstr});

$win = pgwin->new({Device=>'sfdpowerspec.ps/vcps', windowxsize => 8, windowysize => 5});
$win -> line(log($f),$yfit,{Color=> "Red", LineStyle=>2,LineWidth=>2});
$win -> hold;
$win -> line(log($f),log($P),{LineWidth=>2});
$win -> text(' \ga='.sprintf("%6.2f",$index),min(log($f)),average(log($P)));
$win -> release;
#while (!<>) {};





print "hr \n";
&Vtools::view($hr);
print "powerspec \n";
&Vtools::view($powerspec);
$Hr=$hr->copy;
$Hi=$hi->copy;


#$Hr=&restorefft($Hr);
#$Hi=&restorefft($Hi);


#&Vtools::view($Hr);
#$Hr = &restorefft($Hr);
#$Hi = &restorefft($Hi);
#&Vtools::view(exp($Hr/max($Hr)));
$Hmod=sqrt($Hr**2+$Hi**2);
#$Hinv=&restorefft(1/$Hmod);
print "Hmod \n";
&Vtools::view($Hmod);
print "Hr \n";
&Vtools::view($Hr);
print "Hi \n";
&Vtools::view($Hi);
#print "Hinv \n";
#&Vtools::view($Hinv);
#print "exp (Hmod/max) \n";
#&Vtools::view(exp($Hmod/max($Hmod)));



#############################################################
#scaling laws:



$Nbins = 50 ; # $N/2;
$size = zeroes($Nbins);
$Ltot = zeroes($Nbins);
foreach $i (0..$Nbins-1) {
    $size($i) .= ($i+1)*($N/(2*$Nbins))*$dx;
    $Ltot($i) .= sum($Hmod->where($rad<$size(($i))));
}

#&Vtools::spec(log($size),log($Ltot));


$constoff = ones(nelem($size)) ; 
$funcs = cat $constoff,log($size);

$sigs=1/($constoff*0.001*max(log($Ltot)));



($yfit,$coefslin,$covlin) = linfit1d_norm log($Ltot),$funcs,{Weights=>$sigs};

wcols '%5.2e$\pm$%5.2e(%4.1f)', $coefslin,sqrt($covlin->diagonal(0,1)),$coefslin/sqrt($covlin->diagonal(0,1));

$fracdim = sprintf("%5.2f",$coefslin(1) ->sclr);

$win = pgwin->new({Device=>'scale.ps/vps', windowxsize => 8, windowysize => 5});
$win -> line(log($size),$yfit,{Color=> "Red", LineStyle=>2,LineWidth=>2});
$win -> hold;
$win -> line(log($size),log($Ltot),{LineWidth=>2});
$win -> text("  D=$fracdim",min(log($size)),average(log($Ltot)));
$win -> release;
#while (!<>) {};



sub restorefft{
    my $Hr=shift;
    my $N=$Hr->getdim(0);
    $Hr=rotate($Hr,$N/2);
    $Hr=$Hr->xchg(0,1);
    $Hr=rotate($Hr,$N/2);
    $Hr=$Hr->xchg(0,1);

#    rotate($Hr,$N/2);
#    $Hr=$Hr->xchg(0,1);
#    $Hr=rotate($Hr,$N/2);
#    $Hr=$Hr->xchg(0,1);
    
    return $Hr;
}

##p (($fwhm/(2*sqrt(2*log(2))))*2*acos(0)*2*$sigma)