Initial commitHEADmaster

1 files changed, 377 insertions, 0 deletions

diff --git a/harmoneff/harmoneff.pas b/harmoneff/harmoneff.pas
new file mode 100644
index 0000000..20b273a
--- /dev/null
+++ b/harmoneff/harmoneff.pas
@@ -0,0 +1,377 @@
+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.