#!/usr/bin/env perl
# Calculate the self-energy and the spectral function via the self-energy trick
# General approach, for QS symmetry type, for the 2x2 matrix case
# Rok Zitko, rok.zitko@ijs.si, Oct 2008, June 2012, May 2019

use warnings;
use strict;
use Math::Complex;

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

my @imG11 = readcol("c-imG-11.dat", 2);
my @reG11 = readcol("c-reG-11.dat", 2);
my @imF11 = readcol("c-imF-11.dat", 2);
my @reF11 = readcol("c-reF-11.dat", 2);
my @imG22 = readcol("c-imG-22.dat", 2);
my @reG22 = readcol("c-reG-22.dat", 2);
my @imF22 = readcol("c-imF-22.dat", 2);
my @reF22 = readcol("c-reF-22.dat", 2);
my @imG12 = readcol("c-imG-12.dat", 2);
my @reG12 = readcol("c-reG-12.dat", 2);
my @imF12 = readcol("c-imF-12.dat", 2);
my @reF12 = readcol("c-reF-12.dat", 2);
my @imG21 = readcol("c-imG-21.dat", 2);
my @reG21 = readcol("c-reG-21.dat", 2);
my @imF21 = readcol("c-imF-21.dat", 2);
my @reF21 = readcol("c-reF-21.dat", 2);

my @imdelta1 = readcol("Delta-1.dat", 2); # not prefixed!
my @redelta1 = readcol("Delta-1-re.dat", 2);
my @imdelta2 = readcol("Delta-2.dat", 2); # not prefixed!
my @redelta2 = readcol("Delta-2-re.dat", 2);

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

my $imsigma11 = "imsigma-11.dat";
open(IM11, ">$imsigma11") or die "Can't open $imsigma11 for writing, stopped";

my $resigma11 = "resigma-11.dat";
open(RE11, ">$resigma11") or die "Can't open $resigma11 for writing, stopped";

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

my $imsigma22 = "imsigma-22.dat";
open(IM22, ">$imsigma22") or die "Can't open $imsigma22 for writing, stopped";

my $resigma22 = "resigma-22.dat";
open(RE22, ">$resigma22") or die "Can't open $resigma22 for writing, stopped";

# my $pi = pi();

my $len = @omega;
for (my $i = 0; $i < $len; $i++) {
  my $o = $omega[$i];
  
  # Up to a factor of -pi, which however cancels out!
  my $G11 = $reG11[$i] + i * $imG11[$i];
  my $F11 = $reF11[$i] + i * $imF11[$i];
  my $G22 = $reG22[$i] + i * $imG22[$i];
  my $F22 = $reF22[$i] + i * $imF22[$i];
  my $G12 = $reG12[$i] + i * $imG12[$i];
  my $F12 = $reF12[$i] + i * $imF12[$i];
  my $G21 = $reG21[$i] + i * $imG21[$i];
  my $F21 = $reF21[$i] + i * $imF21[$i];
  
  my $det = $G11 * $G22 - $G12 * $G21;
  my $sigma11 = -($F12 * $G21 - $F11 * $G22)/$det;
  my $sigma12 = ($F12 * $G11 - $F11 * $G12)/$det;
  my $sigma21 = -($F22 * $G21 - $F21 * $G22)/$det;
  my $sigma22 = ($F22 * $G11 - $F21 * $G12)/$det;
    
  my $delta1 = i * $imdelta1[$i] + $redelta1[$i];
  my $delta2 = i * $imdelta2[$i] + $redelta2[$i];
    
  my $m11 = $o + $delta1 - $sigma11;
  my $m12 = -$sigma12;
  my $m21 = -$sigma21;
  my $m22 = $o + $delta2 - $sigma22;
    
  my $detm = $m11 * $m22 - $m12 * $m21;
  my $gf11 = $m22/$detm;
  my $gf12 = -$m12/$detm;
  my $gf21 = -$m21/$detm;
  my $gf22 = $m11/$detm;
 
  my $aw11 = -1.0/pi * Im($gf11);
  print SELF11 "$o $aw11\n";
  my $aw22 = -1.0/pi * Im($gf22);
  print SELF22 "$o $aw22\n";

  my $re11 = Re($sigma11);
  my $im11 = Im($sigma11);
  print IM11 "$o $im11\n";
  print RE11 "$o $re11\n";
  my $re22 = Re($sigma22);
  my $im22 = Im($sigma22);
  print IM22 "$o $im22\n";
  print RE22 "$o $re22\n";
}

# Read a column from a file
sub readcol {
    my ($filename, $column) = @_;
    my @l;
    open(F, "<$filename") or die "Can't open $filename for reading, stopped";
    while (<F>) {
	chomp;
	my $x = (split)[$column-1];
	push(@l, $x);
    }
    close(F);
    @l
};