9.9 Unit KoherU

{-------------------------------------------------------------------}

{ Berechnung des Kohaerenzgrades und der Kohaerenzlaenge }
unit KoherU;

{-------------------------------------------------------------------}

interface

procedure gammaBerechnung;
procedure LCBerechnung;

{-------------------------------------------------------------------}

implementation

uses
  Glob,        { Globale Deklarationen }
  LoadSave,    { Laden und Speichern von Dateien }
  TransU,      { Vorbereitungen zur (inversen) FFT }

  { Turbo Pascal Standard Units }
  App, Crt, Dialogs, Objects, Views;

{-------------------------------------------------------------------}

{ Berechnung des Kohaerenzgrades }
procedure gammaBerechnung;

var
  R                                                 : TRect;
  DD                                                : PDialog;
  i, j, Messung, max                                : word;
  I0, I0n, Imax, Imin, Imin1, M, gam, Summe, QSumme : real;
  Iref, Iobj                                        : ^Realarray;
  fnsave, s, s1                                     : string;
  Taste                                             : char;

begin
  New(Iref);
  New(Iobj);
  Taste := ' ';
  fnsave := FileName;
  Path := copy(FileName, 0, Length(FileName) - 4);
  repeat
    Path := copy(Path, 0, Length(Path) - 1);
    s := copy(Path, Length(Path), Length(Path));
  until s = '\';

  for j := 0 to Anzahldelta - 1 do
  begin
    Str(j + 1, s);
    Str(Anzahldelta, s1);
    R.Assign(0, 0, 40, 10);
    DD := New(PDialog, Init(R, 'Berechnung'));
    with DD^ do
    begin
      Options := Options or ofCentered;
      R.Grow(-1, -1);
      Dec(R.B.Y, 3);
      Insert(New(PStaticText, Init(R, #13 + #3 + s + '. von ' + s1
        + ' Abstaenden' + #13 + #13 + #13 + #3
        + 'Abbruch mit <ESC>')));
    end;
    Desktop^.Insert(DD);
    I0 := 0;
    for i := 0 to nmax - 1 do
    begin
      Iref^[i] := 0;
      Iobj^[i] := 0;
    end;

    { Nullintensitaet }
    for Messung := 0 to AnzahlMessung - 1 do
    begin
      FileName := Path + df[j].dark;
      if AnzahlMessung > 1 then
      begin
        Str(Messung, s);
        FileName := FileName + '_' + s;
      end;
      FileName := FileName + '.' + ExtMess;
      LoadFileI(Imes^, FileName);
      I0n := 0;
      for i := 0 to nmax - 1 do I0n := I0n + Imes^[i];
      I0n := I0n / nmax;
      I0 := (I0n + Messung * I0) / (Messung + 1);
    end;

    { Referenzstrahl }
    for Messung := 0 to AnzahlMessung - 1 do
    begin
      FileName := Path + df[j].ref;
      if AnzahlMessung > 1 then
      begin
        FileName := FileName + '_' + s;
        Str(Messung, s);
      end;
      FileName := FileName + '.' + ExtPlot;
      Synthese;
      for i := 0 to nmax - 1 do
      begin
        Itrans^[i] := Itrans^[i] - I0;
        Iref^[i] :=
          Iref^[i] + (Itrans^[i] - Iref^[i]) / (Messung + 1);
      end;
    end;

    { Objektstrahl }
    for Messung := 0 to AnzahlMessung - 1 do
    begin
      FileName := Path + df[j].obj;
      if AnzahlMessung > 1 then
      begin
        FileName := FileName + '_' + s;
        Str(Messung, s);
      end;
      FileName := FileName + '.' + ExtPlot;
      Synthese;
      for i := 0 to nmax - 1 do
      begin
        Itrans^[i] := Itrans^[i] - I0;
        Iobj^[i] :=
          Iobj^[i] + (Itrans^[i] - Iobj^[i]) / (Messung + 1);
      end;
    end;

    { Interferenzmuster }
    w[j].gamma := 0;
    Summe := 0;
    QSumme := 0;
    for Messung := 0 to AnzahlMessung - 1 do
    begin
      FileName := Path + df[j].int;
      if AnzahlMessung > 1 then
      begin
        FileName := FileName + '_' + s;
        Str(Messung, s);
      end;
      FileName := FileName + '.' + ExtPlot;
      Synthese;

      { Maximalintensitaet }
      Imax := 0;
      for i := 0 to nmax - 1 do
        if Itrans^[i] > Imax then
        begin
          Imax := Itrans^[i];
          max := i;
        end;

      { Minimalintensitaet }
      Imin := Imax;
      if (max > 0.8 * Tau) and (max < nmax - 1 - 0.8 * Tau) then
      begin
        for i := max downto max - round(0.75 * Tau) do
          if Itrans^[i] < Imin then Imin := Itrans^[i];
        Imin1 := Imax;
        for i := max to max + round(0.75 * Tau) do
          if Itrans^[i] < Imin1 then Imin1 := Itrans^[i];
        Imin := (Imin + Imin1) / 2;
      end
      else
      begin
        if max > 0.8 * Tau then
        begin
          for i := max downto max - round(0.75 * Tau) do
            if Itrans^[i] < Imin then Imin := Itrans^[i];
        end
        else
          for i := max to max + round(0.75 * Tau) do
            if Itrans^[i] < Imin then Imin := Itrans^[i];
      end;

      Imax := Imax - I0;
      Imin := Imin - i0;

      { Modulation }
      M := (Imax - Imin) / (Imax + Imin);

      { Kohaerenzgrad }
      gam := M * (Iref^[max] + Iobj^[max]) /
        (2 * sqrt(Iref^[max] * Iobj^[max]));

      { Mittelwert ueber die Messungen }
      w[j].gamma :=
        (gam + Messung * w[j].gamma) / (Messung + 1);

      { Standardabweichung }
      Summe := Summe + gam;
      QSumme := QSumme + sqr(gam);
      if Messung > 0 then
        w[j].SAM := sqrt(abs((QSumme - 2 * Summe * w[j].gamma
          + (Messung + 1) * sqr(w[j].gamma)) / (Messung + 1)
          * Messung))
      else w[j].SAM := 0;
    end;
    if keypressed then
    begin
      Taste := readkey;
      if Taste = chr(27) then Abbruch := true;
    end;
    Dispose(DD, Done);
    if Abbruch then j := Anzahldelta - 1;
  end;
  FileName := fnsave;
  Dispose(Iref);
  Dispose(Iobj);
end; { gammaBerechnung }

{-------------------------------------------------------------------}

{ Berechnung der Kohaerenzlaenge }
procedure LCBerechnung;

var
  i                         : word;
  dc, gam, Summe, tmp, tmpa : real;

begin

  { Kohaerenzlaenge }
  tmp := 1;
  for i := 0 to Anzahldelta - 1 do
    if abs(w[i].gamma - exp(0 - 1)) < tmp then
    begin
      deltac := w[i].delta;
      tmp := abs(w[i].gamma - exp(0 - 1));
    end;

  Summe := 0;
  for i := 0 to Anzahldelta - 1 do
  begin
    gam := exp(0 - sqr(w[i].delta / deltac));
    Summe := Summe + gam - w[i].gamma;
  end;
  tmp := Summe / Anzahldelta;
  dc := deltac / 1000;

  repeat
    deltac := deltac + dc;
    tmpa := tmp;
    Summe := 0;
    for i := 0 to Anzahldelta - 1 do
    begin
      gam := exp(0 - sqr(w[i].delta / deltac));
      Summe := Summe + gam - w[i].gamma;
    end;
    tmp := Summe / Anzahldelta;
  until abs(tmp) > abs(tmpa);
  deltac := deltac - dc;

  repeat
    deltac := deltac - dc;
    tmpa := tmp;
    Summe := 0;
    for i := 0 to Anzahldelta - 1 do
    begin
      gam := exp(0 - sqr(w[i].delta / deltac));
      Summe := Summe + gam - w[i].gamma;
    end;
    tmp := Summe / Anzahldelta;
  until abs(tmp) > abs(tmpa);
  deltac := deltac + dc;

  { Standardabweichung }
  Summe := 0;
  for i := 0 to Anzahldelta - 1 do
  begin
    tmp := w[i].SAM / w[i].gamma;
    Summe := Summe + tmp;
  end;
  deltacSAM := deltac * Summe / Anzahldelta;
end; { LCBerechnung }

{-------------------------------------------------------------------}

end. { KoherU }

{-------------------------------------------------------------------}