#!/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);

$N=512;
$dx=(1/($N-1));
$h=zeroes($N,$N);
$x=sequence($N)*$dx;
$xs=$dx*($h->xvals);
$ys=$dx*($h->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);
$Hausdorf_1D  = 0.8  ;
$index  = 1*(-1-2*$Hausdorf_1D);
$index  = -3.7;
print "index is $index \n";
#print "expect a 2d fractal dim of ",2*$Hausdorf_1D," \n";
$k0=1;
#$hr=
#$hi=(random($N,$N)-0.5)*((($rad+$dx)/$k0)**($index/2));



$amp = (grandom($N,$N)-0.5)*((($rad+$dx)/$k0)**($index/2));

#$amp($N/2,$N/2) = 0;

$phase = 2*$pi*(random($N,$N));

$subim = $amp($N/2:$N-1,$N/2:$N-1);

#$tr = t_linear({rot=>180,origin=>[$N/2-1,$N/2-1]});
$tr = t_linear({rot=>180});
#&Vtools::view($subim);
#$subim2 = $subim -> map ($tr);
#&Vtools::view($subim2);

$symamp = $amp -> copy;
$symphase = $phase -> copy;
$ramp = $amp -> map($tr) -> copy;
$rphase = $phase -> map($tr) -> copy;

&Vtools::view($amp);


###symmetrisation about one pixel
$symamp(1:$N/2,1:$N/2) .= $amp($N/2:$N-1,$N/2:$N-1) -> map($tr);
$symphase(1:$N/2,1:$N/2) .= -1 * $phase($N/2:$N-1,$N/2:$N-1) -> map($tr);
$symamp($N/2+1:$N-1,1:$N/2-1) .= $amp(1:$N/2-1,$N/2+1:$N-1) -> map($tr);
$symphase($N/2+1:$N-1,1:$N/2-1) .= -1 * $phase($N/2+1:$N-1,1:$N/2-1) -> map($tr);
#

#
####symmetrisation about one intersection
#$symamp(0:$N/2-1,0:$N/2-1) .= $amp($N/2:$N-1,$N/2:$N-1) -> map($tr);
#$symphase(0:$N/2-1,0:$N/2-1) .= -1 * $phase($N/2:$N-1,$N/2:$N-1) -> map($tr);
#
#$symamp($N/2:$N-1,0:$N/2-1) .= $amp(0:$N/2-1,$N/2:$N-1) -> map($tr);
#$symphase($N/2:$N-1,0:$N/2-1) .= -1 * $phase(0:$N/2-1,$N/2:$N-1) -> map($tr);
#
##old> #
##old> #$symamp(0:$N/2-1,0:$N/2-1) .= $ramp(0:$N/2-1,0:$N/2-1); 
##old> #$symphase(0:$N/2-1,0:$N/2-1) .= -1 * $rphase(0:$N/2-1,0:$N/2-1)  ;
##old> #$symamp($N/2:$N-1,0:$N/2-1) .= $ramp($N/2:$N-1,0:$N/2-1);
##old> #$symphase($N/2:$N-1,0:$N/2-1) .= -1 * $rphase ($N/2:$N-1,0:$N/2-1);
##old> #
##old> 
##old> #$symamp(1:$N/2,1:$N/2) .= $ramp(0:$N/2-1,0:$N/2-1); 
##old> #$symphase(1:$N/2,1:$N/2) .= -1 * $rphase(0:$N/2-1,0:$N/2-1)  ;
##old> #$symamp($N/2:$N-1,1:$N/2) .= $ramp($N/2:$N-1,0:$N/2-1);
##old> #$symphase($N/2:$N-1,1:$N/2) .= -1 * $rphase ($N/2:$N-1,0:$N/2-1);
##old> #

$amp = $symamp;
$phase = $symphase;
#
print "SYMMETRY TEST:  \n";

#print "10,10 :",$phase($N/2+9,$N/2+9),"\n";
#print "10,10 :",$phase($N/2-10,$N/2-10),"\n";
#print "10,10 :",$amp($N/2+9,$N/2+9),"\n";
#print "10,10 :",$amp($N/2-10,$N/2-10),"\n";



#&Vtools::view($amp);
#&Vtools::view($phase);


$hr=$amp*cos($phase);
$hi=$amp*sin($phase);

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

$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)+($N/(2*$Nbins))*$dx) )));
}
$xtext=max($f)/2.;
$ytext=max($P)/2.;
$codstr ='$win->text(\' \ga='.$index."\',$xtext,$ytext,{CHARSIZE=>2})";
print "CODE $codstr \n";
&Vtools::spec($f,$P,{code=>$codstr});


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

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

ifftnd($Hr,$Hi);

#$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)

