// Long's correction of Raman spectra and normalisation
// last rev. Oct 2010
// ensures strictly increasing values of wavenumber
// calc. e.s.e. as Long cor. norm. sqrt(n_raw) 3d output col.
// exp. format to avoid null ese.
program long3;

uses getopts, mtx3;

type
  func_dd = function(x:double):double;

var
  MINESE: double = 1.0; // KLUDGE FOR BAD DATA with 0 or <0 counts !
  // for normal data ese is sqrt(counts). It happens that 0 or negative
  // counts are reported by the brain damaged software; so I enforce
  // as a workaround, that the minimum ese must be 1.0 so that these
  // data can be processed by gcv spline smoothing. Option -e can
  // overwrite this value
  nu0: double = 1.943513717e6; // m^(-1) wave number  lambda=514.532 nm
  T: double = 273.15+25.0;     // K  default temperature in Kelvin


// Note that system.sqrt prototype is special and must be wrapped.
// prototype of sqrt(): function(Extended):Extended;InternProc
// It cannot be used for custom functions... AFAIK.
function dsqrt(x:double):double;
begin
  if x>0 then
    dsqrt := sqrt(x)
  else
    dsqrt := MINESE // BAD KLUDGE FOR BRAIN DAMAGE DATA
end;

// this belongs to mtx3...
procedure fct(var r:mtx; a:mtx; f:func_dd);
var
  i, j: longint;
begin
  zer(r,lgn(a),col(a));
  for i := 1 to lgn(a) do begin
    for j := 1 to col(a) do begin
      elm(r,i,j,f(elm(a,i,j)))
    end
  end
end;

procedure help;
begin
  writeln;
  writeln('Raman spectra correction and normalisation using Long''s formula.');
  writeln('Assuming e.s.e. on raw counts to be sqrt(counts) then corrected e.s.e. are returned');
  writeln;
  writeln('roux@ipgp.jussieu.fr');
  writeln;
  writeln('USAGE: long3 < raw.xy > lgn.xy');
  writeln;
  writeln('OPTIONS');
  writeln;
  writeln('  -l wavelength in nm; def: 514.532');
  writeln('  -t temperature in Celsius; def: 25.0');
  writeln('  -e KLUDGE e.s.e. for counts <= 0; def: 1.0');
  writeln;
  halt
end;

procedure options;
var
  code: integer;
  c: char;
  v: double;
  optstr: string = 'hl:t:e:';
begin
  c := getopt(optstr);
  while c<>EndOfOptions do begin
    case c of
      'h': help;
      'l': begin // lambda en nm
             val(optarg,v,code);
             assert(code=0);
             nu0 := 1.0/v*1e9;
           end;
      't': begin
             val(optarg,v,code);
             assert(code=0);
             T := 273.15 + v
           end;
      'e': begin
             val(optarg,v,code);
             assert(code=0);
             MINESE := v
           end;
    end;
    c := getopt(optstr)
  end
end;

function long(x, y:double):double;
const
  h: double = 6.62606896e-34;   // J.s    Plank
  k: double = 1.38066e-23;      // J/K    Boltzman
  c: double = 2.9979e8;         // m/s    Speed of light
var
  rnu: double; // nu0-nu        m^(-1)
  nu: double;  // nu            m^(-1)
  t0, t1, t2: double;
begin
  nu := 100.0*x; // cm-1 -> m-1 Raman shift
  rnu := nu0-nu; // nu0 is in m-1
  t0 := nu0*nu0*nu0*nu/rnu/rnu/rnu/rnu;
  t1 := -h*c*nu/k/T; // c in m/s  : t1 dimensionless
  t2 := 1 - exp(t1);
  long := double(y*t0*t2)
end;

procedure long(var r: mtx; x, y: mtx);
var
  xi, yi, zi: double;
  i, n: longint;
begin
  n := lgn(x);
  assert( lgn(y)=n );
  assert( col(x)=1 );
  assert( col(y)=1 );
  zer(r,n,1);
  for i := 1 to n do begin
    xi := elm(x,i,1);
    yi := elm(y,i,1);
    zi := long(xi,yi);
    elm(r,i,1,zi)
  end
end;

var
  d: mtx;
  x, y, w, e, z: mtx;
  i: longint;
  x0, x1, y0, y1, e0, e1, norm: double;

begin
  options;
  mrd(d);
  x := col(d,1);
  y := col(d,2);
  assert( min(y)>0 );
  fct(w,y,@dsqrt); // w = sqrt(y)
  long(z,x,y);
  norm := max(z);
  long(e,x,w);
  dvd(z,norm); // z = z / norm
  dvd(e,norm); // e = e / norm
// ensure precision and strict monotony of x
// this is why I don't use mwr() for output
  x0 := elm(x,1,1);
  y0 := elm(z,1,1);
  e0 := elm(e,1,1);
  writeln(x0:12:3, y0:12:6, e0:12);
  for i := 2 to lgn(d) do begin
    x1 := elm(x,i,1);
    y1 := elm(z,i,1);
    e1 := elm(e,i,1);
    if ( (x1-x0) > 0.1 ) then begin
      writeln(x1:12:3, y1:12:5, e1:12)
    end;
    x0 := x1
  end
end.
