#!/usr/bin/env perl
# Calculate the spectral function via the self-energy trick
# Rok Zitko, zitko@theorie.physik.uni-goettingen.de, Oct 2008

use warnings;
use Math::Complex;
require "misc.pl";

my $U = `getparam U param.loop`; # band e-e repulsion
my $delta = `cat param.delta`;
my $ed = $delta - $U/2;

chomp($U, $delta);
print "U=$U delta=$delta ed=$ed\n";

# Hybridisation function
system "kk Delta.dat Delta-re.dat";

chdir "res";

my @omega = readcol("c-imG.dat", 1);

my @imG = readcol("c-imG.dat", 2);
my @reG = readcol("c-reG.dat", 2);
my @imF = readcol("c-imF.dat", 2);
my @reF = readcol("c-reF.dat", 2);

my @imdelta = readcol("../Delta.dat", 2);
my @redelta = readcol("../Delta-re.dat", 2);

my $selffn = "c-self.dat";
open(SELF, ">$selffn") or die "Can't open $selffn for writing, stopped";

my $len = @omega;
for ($i = 0; $i < $len; $i++) {
  my $o = $omega[$i];
  my $G = $reG[$i] + i * $imG[$i];

  my $F = $reF[$i] + i * $imF[$i];
  my $sigmaU = $U * $F/$G;
    
  my $sigma = $sigmaU;
  my $delta = i * $imdelta[$i] + $redelta[$i];
 
  my $gf = 1.0/($o - $ed + $delta - $sigma);
  my $aw = -1.0/pi * Im($gf);

  print SELF "$o $aw\n";
}
