path: root/harmoneff/harmoneff.pas

program harmoneff;

uses sysutils, myMath, classes, process, math;

Type TTodo = record
              p_pol,
              transmitted,
              plasma:      boolean;
              order:       integer;
             end;
     TLookAt = record
                key,name: string;
               end;

var datadir: string;

function catchValue(fnam,key,name: string): string;
var f:  textfile;
    s:  string;
    kf: boolean;
begin
 if not fileexists(fnam) then Writeln(fnam);
 assign(f,fnam);
 reset(f);
 kf:=false;
 while not eof(f) do
  begin
   readln(f,s);
   if pos('#',s)>0 then s:= copy(s,1,pos('#',s)-1);
   while pos(' ',s)>0 do
    delete(s,pos(' ',s),1);
   while pos(#9,s)>0 do
    delete(s,pos(#9,s),1);
   if length(s)=0 then continue;
   if s[1]='&' then
    begin
     kf:='&'+key=s;
     continue;
    end;
   if kf and (pos(name+'=',s)=1) then
    begin
     delete(s,1,length(name)+1);
     catchValue:=s;
     close(f);
     exit;
    end;
  end;
 close(f);
 catchValue:='-1';//'arrrgh, couldn''t find &'+key+' -> '+name+' in '+fnam+'!';
end;

function getIntensity(wen: string; Spalte: Integer; Ordnung,Otol: extended; fit: boolean): TExtPoint;
var F:         Textfile;
    S:         String;
    n:         float;
    I,J:       integer;
    Daten:     array of TExtPoint;
    y0,dt,err: extended;
    erg:       TExtPoint;
begin
 if not fileexists(datadir+'/Post/ft-'+wen) then writeln(datadir+'/Post/ft-'+wen);
 assign(F,datadir+'/Post/ft-'+wen);
 reset(F);
 n:=0;
 Setlength(Daten,0);
 while not eof(F) do
  begin
   readln(F,S);
   while (length(S)>0) and (S[1] in [#9,' ']) do
    delete(S,1,1);
   if (length(S)=0) or (S[1]='#') then continue;
   S:=S+' ';
   n:=strtofloat(copy(S,1,pos(' ',S)-1));
   for I:=1 to Spalte do
    begin
     delete(S,1,pos(' ',S));
     while (length(S)>0) and (S[1] in [#9,' ']) do
      delete(S,1,1);
    end;
   S:=copy(S,1,pos(' ',S)-1);
   if (abs(n-Ordnung)<=Otol/2) or (Otol<0) then
    begin
     Setlength(Daten,length(Daten)+1);
     Daten[length(Daten)-1].x:=n-Ordnung;
     Daten[length(Daten)-1].y:=strtofloat(S);
    end;
  end;
 close(F);
 err:=0.001;

 if length(Daten)=0 then writeln('keine Daten zum fitten!');

 if fit then
  begin
   y0:=0; // generate some usable values of y0 and dt to start with
   I:=0;
   while (Daten[I].x<0) and (I<length(Daten)) do
    begin
     y0:=Daten[I].y;
     inc(I);
    end;
   J:=I;
   while (I<length(Daten)-1) and (Daten[I].y > exp(-1)*y0) do
    inc(I);
   while (J>0) and (Daten[J].y > exp(-1)*y0) do
    dec(J);
   dt:=(Daten[I].x-Daten[J].x)/2;
   
   if y0=0 then writeln('y0 ist vorher null!');
   
   GaussFitDownhillSimplex(Daten,y0,dt,err);
   if y0=0 then writeln('y0 ist nachher null!');
   erg.val:=abs(y0*dt*sqrt(pi));
   erg.err:=abs(sqrt(err/(length(Daten)-1))/y0);
  end
   else
  begin
   erg.val:=0;
   erg.err:=0;
   J:=0;
   for I:=0 to length(Daten)-1 do
    begin
     erg.val:=erg.val+Daten[I].y;
     if (10*I<length(Daten)) or 
        ((length(Daten)-I-1)*10<length(Daten)) then
      begin
       erg.err:=erg.err+Daten[I].y;
       inc(J);
      end;
    end;
   erg.val:=erg.val*(Daten[1].x-Daten[0].x);
   if J>0 then erg.err:=erg.err/J*(Daten[length(Daten)-1].x-Daten[0].x);
   erg.err:=erg.err/erg.val;
  end;
 getIntensity:=erg;
end;

procedure fehler(errnum: integer);
var i: integer;
begin
 case errnum of
  0:
   begin
    writeln('number of given Params was '+inttostr(paramcount)+':');
    for i:=0 to paramcount do
     writeln(' '+inttostr(i)+': "'+paramstr(i)+'"');
    writeln('usage:');
    writeln(' '+paramstr(0)+' output input ($key $name)* \# (p|s)(t|r)(P?[0..9]*)');
    writeln('   -- efficiency of p/s-pol., transm./refl., n-th (plasma) harmonic over $key-$name''s');
    writeln(' '+paramstr(0)+' output ($key $name)* \# (p|s)(t|r)(P?[0..9]*)');
    writeln('   -- headlines for above');
   end;
  1:
   begin
    writeln('error: corrupt harmoneff.ini');
   end;
  else
   writeln('*** unknown error: error '+inttostr(errnum)+' ***');
 end{of Case};
 halt;
end;

function  readString(Key,Nam: string): String;
var f:        textfile;
    s:        string;
    rightKey: boolean;
begin
 assign(f,extractfilepath(Paramstr(0))+'harmoneff.ini');
 reset(f);
 rightKey:=false;
 while not eof(f) do
  begin
   readln(f,s);
   if pos('#',s)>0 then delete(s,pos('#',s),length(s));
   while (length(s)>0) and (s[1] in [' ',#9]) do
    delete(s,1,1);
   while (length(s)>0) and (s[length(s)] in [' ',#9]) do
    delete(s,length(s),1);
   if length(s)=0 then continue;
   if s[1]='&' then
    begin
     delete(s,1,1);
     rightKey:=copy(s,1,length(Key))=Key;
     continue;
    end;
   if not rightKey then continue;
   if copy(s,1,length(Nam))=Nam then
    begin
     delete(s,1,length(Nam));
     delete(s,1,pos('=',s));
     while (length(s)>0) and (s[1] in [' ',#9]) do
      delete(s,1,1);
     readString:=s;
     close(f);
     exit;
    end;
  end;
 close(f);
 writeln('&'+Key+' '+Nam+' not found!');
 fehler(1);
end;

function  readFloat(Key,Nam: string): extended;
begin
 readFloat:=strtofloat(readString(Key,Nam));
end;

function  readInt(Key,Nam: string): integer;
begin
 readInt:=strtoint(readString(Key,Nam));
end;

function  readBool(Key,Nam: string): boolean;
begin
 readBool:=readInt(Key,Nam)<>0;
end;

var inputintens,
    outputintens,
    zuPapier:     TExtPoint;
    i,j,numl:     integer;
    todos:        array of TTodo;
    lookats:      array of TLookAt;
    heading:      boolean;
    S:            string;
    nion,max_err,
    lfbreite,
    pfbreite:     extended;
    FOut:         Textfile;
    fitten,
    mit_fehler,
    plusminus:	  boolean;
    spalte:       array[boolean] of Integer; //  [transm?]
    input_pol:    boolean; // p-pol?
const
    wer:    array[boolean,boolean] of String  = (('gp','fm'),('gm','fp'));  //  [transm?,p-pol?]

procedure readIni;
begin
 fitten:=readBool('integration','fitten');
 mit_fehler:=readBool('output','mit_fehler');
 plusminus:=readBool('output','plusminus');
 max_err:=readFloat('output','max_err');
 spalte[false]:=readInt('daten','l_pos');
 spalte[true]:=readInt('daten','l_pos');
 input_pol:=readBool('daten','input_pol');
 lfbreite:=readFloat('integration','lfbreite');
 pfbreite:=readFloat('integration','pfbreite');
end;

begin
 if paramcount <= 1 then
  fehler(0);
 readIni;
 Assign(FOut,Paramstr(1));
 if fileexists(Paramstr(1)) then
   Append(FOut)
  else
   Rewrite(FOut);
 heading:=not fileexists(Paramstr(2));
 setlength(lookats,0);
 numl:=-1;
 For i:=2+byte(not heading) to Paramcount do
  if (Paramstr(i)='#') or (Paramstr(i)='\#') then 
   begin
    if numl<>-1 then fehler(0);
    numl:=i-1;
   end;
 if (numl=-1) or (odd(numl) xor heading) then fehler(0);
 for i:=1 to (numl-1-Byte(not heading)) div 2 do
  begin
   setlength(lookats,length(lookats)+1);
   lookats[length(lookats)-1].key:=Paramstr(Byte(not heading)+2*i);
   lookats[length(lookats)-1].name:=Paramstr(Byte(not heading)+2*i+1);
  end;
 inc(numl,2);
 
 setlength(todos,0);
 For i:=numl to Paramcount do
  begin
   S:=paramstr(i);
   if (length(S)<3) or
      not (upcase(S[1]) in ['P','S']) or
      not (upcase(S[2]) in ['T','R']) or
      not (upcase(S[3]) in ['P','0'..'9']) then fehler(0);
   for j:=4 to length(S) do
    if not (S[j] in ['0'..'9']) then fehler(0);
   setlength(todos,length(todos)+1);
   todos[length(todos)-1].p_pol:=upcase(S[1]) = 'P';
   todos[length(todos)-1].transmitted:=upcase(S[2]) = 'T';
   todos[length(todos)-1].plasma:=upcase(S[3]) = 'P';
   delete(S,1,2+Byte(todos[length(todos)-1].plasma));
   todos[length(todos)-1].order:=strtoint(S);
  end;
 
 if heading then
  begin
   write(FOut,'source');
   for i:=0 to length(lookats)-1 do
    write(FOut,#9+lookats[i].key+'-'+lookats[i].name);
   for i:=0 to length(todos)-1 do
    begin
     write(FOut,#9);
     if todos[i].p_pol then S:='p'
                       else S:='s';
     S:=S+'-pol ';
     if todos[i].transmitted then S:=S+'transm.'
                             else S:=S+'refl.';
     S:=S+' '+inttostr(todos[i].order);
     if todos[i].plasma then S:=S+' (pl.)';
     write(FOut,S);
     if mit_fehler then write(FOut,#9'd_abs'#9'd_rel');
     if plusminus then write(FOut,#9'min'+S+#9'max'+S);
    end;
   writeln(FOut);
  end
   else
  begin
   datadir:=catchValue(Paramstr(2),'output','path');
   if not directoryexists(datadir) then
    begin
     delete(datadir,1,pos('/',datadir));
     if directoryexists('../daten_alt/'+datadir) then
      datadir:='../daten_alt/'+datadir
       else
      begin
       if directoryexists('../daten_mittelalt/'+datadir) then
        datadir:='../daten_mittelalt/'+datadir
         else
        begin
         writeln('**************************************************************************');
         writeln(datadir);
         halt;
        end;
      end;
    end;
   inputintens:=getIntensity(wer[true,input_pol],1,1,-1,fitten); //'fp'
   write(FOut,datadir);
   for i:=0 to length(lookats)-1 do
    write(FOut,#9+myfloattostr(strtofloat(catchValue(Paramstr(2),lookats[i].key,lookats[i].name))));
   for i:=0 to length(todos)-1 do
    begin
     if todos[i].plasma then
      begin
       nion:=strtofloat(catchValue(Paramstr(2),'box','n_ion_over_nc'));
       outputintens:=getIntensity(wer[todos[i].transmitted,todos[i].p_pol],
                                  Spalte[todos[i].transmitted],
                                  todos[i].order*sqrt(nion),
                                  pfbreite/2*sqrt(nion),fitten);
      end
       else
      begin
       outputintens:=getIntensity(wer[todos[i].transmitted,todos[i].p_pol],
                                  Spalte[todos[i].transmitted],
                                  todos[i].order,
                                  lfbreite,fitten);
      end;
     zuPapier:=divide(outputintens,inputintens);
     if (zuPapier.Err <= max_err) or (max_err < 0) then
      begin 
       write(FOut,#9+valuetostr(zuPapier,mit_fehler));
       if plusminus then
        write(FOut,#9+myfloattostr(zuPapier.Val*(1-zuPapier.Err))+
                   #9+myfloattostr(zuPapier.Val*(1+zuPapier.Err)));
      end
       else
      begin 
       for J:=0 to 2*Byte(mit_fehler) + 2*Byte(plusminus) do
        write(FOut,#9);
      end;
    end;
   writeln(FOut);
  end;

 Close(FOut);

end.