#!/bin/sh
math --batchinput <<END

gap = 0.01;

do[sign_: 1] := Module[{},
   l = BinaryReadList["spec_CFS_dens_A_d-A_d.bin", 
     "Real64"];
   If[l == \$Failed, Exit[]];
   
   l = Partition[l, 2];

   (* Positive or negative part of the spectrum! *)
   lp = Select[l, #[[1]]*sign > 0 &];

   (* Subgap part of the spectrum *)
   l2 = Select[lp, Abs[#[[1]]] < gap &];

   (* Total subgap spectral weight *)
   tot2 = Total[l2[[All, 2]]];
   Print["total subgap weight=", tot2];

   (* Accurate energy and weight from delta peak series *)
   enwe[l_] := Module[{},
     totweight = Total[l[[All, 2]]];
     avgenergy = l[[All, 1]].l[[All, 2]]/totweight;
     {avgenergy, totweight}
     ];

   (* Keep large weight parts, i.e. subgap excitations *)
   l3 = Select[l2, #[[2]]/tot2 > 10^-6 &];

   (* Length along the energy axis *)
   df[{e1_, w1_}, {e2_, w2_}] := Abs[e1 - e2];

   (* Cluster delta peaks according to their energy separation. 
   This allows to separate contributions to different excitations. *)
   cl = FindClusters[l3, DistanceFunction -> df];

   (* Calculate energy means and total weights *)
   result = Map[enwe, cl];

   (* Rescale in units of the gap *)
   result[[All, 1]] = result[[All, 1]]/gap;
   
   result = Select[result, #[[2]] > 10^-3 tot2 && #[[2]] > 10^-6 &];

   Print["result=", result];
   result = Append[result, {}];
   Print[Export[ToString[sign] <> ".dat", result]];
];

do[1];
do[-1];
END
