{-----------------------------------------------------------------------------
 Unit Name: Main
 Author:    kaetzel
 Purpose:   Contains all event handling procedures of the main window
-----------------------------------------------------------------------------}

unit Main;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls,JvDSADialogs, ExtCtrls, JvEdit, Grids,JclMath,JclBase,
  ComplexNumbers,ComplexLinSysSolve, JvCombobox;

type
  TMainWindow = class(TForm)
    AggrConfigButton:   TButton;
    CloseButton:        TButton;
    ResetButton:        TButton;
    Label1:             TLabel;
    Panel1:             TPanel;
    Label2:             TLabel;
    Panel2:             TPanel;
    Label3:             TLabel;
    Label4:             TLabel;
    RefrIndexInput:     TJvEdit;
    WavelengthInput:    TJvEdit;
    StartButton:        TButton;
    Panel3:             TPanel;
    Panel4:             TPanel;
    Panel5:             TPanel;
    Label5:             TLabel;
    Label6:             TLabel;
    Label12:            TLabel;
    Label13:            TLabel;
    EdExtCoeff:         TJvEdit;
    EdScatCoeff:        TJvEdit;
    EdAbsCoeff:         TJvEdit;
    EdAsyPar:           TJvEdit;
    Label14:            TLabel;
    Label15:            TLabel;
    Label16:            TLabel;
    Label17:            TLabel;
    Label18:            TLabel;
    Label19:            TLabel;
    Label20:            TLabel;
    Label21:            TLabel;
    Label22:            TLabel;
    Label23:            TLabel;
    Label24:            TLabel;
    Label25:            TLabel;
    Label26:            TLabel;
    Label27:            TLabel;
    Label28:            TLabel;
    Label29:            TLabel;
    Label30:            TLabel;
    ApAngleInput:       TJvEdit;
    Label31:            TLabel;
    Label32:            TLabel;
    Label33:            TLabel;
    Label34:            TLabel;
    Label35:            TLabel;
    Label36:            TLabel;
    Label37:            TLabel;
    Label38:            TLabel;
    Label39:            TLabel;
    Label40:            TLabel;
    CbOrientations:     TComboBox;
    Label41:            TLabel;
    Label7: TLabel;
    Label8: TLabel;
    Label9: TLabel;
    Label10: TLabel;
    Label11: TLabel;
    Label42: TLabel;
    EdRadPresCoeff: TJvEdit;
    Label43: TLabel;
    Label44: TLabel;
    Label45: TLabel;
    Label46: TLabel;
    procedure AggrConfigButtonClick(Sender: TObject);
    procedure CloseButtonClick(Sender: TObject);
    procedure ResetButtonClick(Sender: TObject);
    procedure StartButtonClick(Sender: TObject);
    procedure FormDestroy(Sender: TObject);
    procedure FormCreate(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var                               //here all global variables that are used in
  MainWindow:         TMainWindow;//nearly all units of the program are declared
  ButtonClicked:      bool;              //Inquiry whether Button is clicked
  CloseClicked:       bool;
  Nl:                 Word;              //Number of Spheres in agglomerate
  dia:                array of real;     //diameter of the spheres
  x:                  array of extended; //size parameter of the spheres
  Xcoord:             array of float;    //cartesian coordinates of the spheres
  Ycoord:             array of float;
  Zcoord:             array of float;
  Rcoord:             array of float;    //Spherical coordinates of the spheres
  Thetacoord:         array of float;
  Phicoord:           array of float;
  RefrIndexReal:      array of real;     //Refractive Indices of the spheres
  RefrIndexImag:      array of real;
  RefrIndMed,WLength: real;              //Refractive Index of the medium and wavelength
  WiscombeNumber:     array of SmallInt; //Largest Wiscombe Criterion of the spheres
  RelRefrIndex:       array of array of real; //Relative refractive index
  MieParameter:       array of array of array of extended; //a(n) and b(n) for each sphere
  Wavenumber:         float;             // 2*Pi/lambda
  djl:                array of array of float; //Distance Matrix coordinates
  phijl:              array of array of float;
  thetajl:            array of array of float;
  Vector_pq:          TComplexVector;  //Vector of primary expansion coeff's
  SolutionVector:     TComplexVector;
  Matrix_AB:          TComplexMatrix;  //Vector translation coefficient matrix
  Vector_amn_bmn:     TComplexVector;  //Vector of partial scattering coefficients
  Vector_amnl_bmnl:   TComplexVector;
  X0coord:            array of float;    //Coordinates of first orientation
  Y0coord:            array of float;
  Z0coord:            array of float;

implementation

uses NumSpheres, Coordinates, Functions, MieParam;

{$R *.dfm}


{-----------------------------------------------------------------------------
  Procedure: TMainWindow.AggrConfigButtonClick
  Purpose:   Calls the agglomerate configuration windows
  Comments:
-----------------------------------------------------------------------------}

procedure TMainWindow.AggrConfigButtonClick(Sender: TObject);
begin
  AggrConfigButton.Enabled:=False;
        //"Agglomerate configuration.."
  Label22.Show;
  NumSpheresWindow.Show;
end;


{-----------------------------------------------------------------------------
  Procedure: TMainWindow.CloseButtonClick
  Purpose:   Closes the application
  Comments:
-----------------------------------------------------------------------------}

procedure TMainWindow.CloseButtonClick(Sender: TObject);
begin
  CloseClicked:=True;
  Close;
end;


{-----------------------------------------------------------------------------
  Procedure: TMainWindow.ResetButtonClick
  Purpose:   Reinitialization of program
  Comments:
-----------------------------------------------------------------------------}

procedure TMainWindow.ResetButtonClick(Sender: TObject);
var
  i,j:integer;
  buttons:array[1..2] of string;
  results:array[0..1] of integer;
  ResultHandler:TModalResult;
begin
  SetLength(MieParameter,0,0,0);
  SetLength(Vector_pq,0);
  SetLength(Matrix_AB,0,0);
  SetLength(Vector_amn_bmn,0);
  SetLength(Vector_amnl_bmnl,0);
  buttons[1]:='New Agglomerate';
  buttons[2]:='Other Parameters';
  results[0]:=1;
  results[1]:=2;
  EdExtCoeff.Enabled:=True;
  EdScatCoeff.Enabled:=True;
  EdAbsCoeff.Enabled:=True;
  EdAsyPar.Enabled:=True;
  EdRadPresCoeff.Enabled:=True;
  EdExtCoeff.Text:='';
  EdScatCoeff.Text:='';
  EdAbsCoeff.Text:='';
  EdAsyPar.Text:='';
  EdRadPresCoeff.Text:='';
  EdExtCoeff.Enabled:=False;
  EdScatCoeff.Enabled:=False;
  EdAbsCoeff.Enabled:=False;
  EdAsyPar.Enabled:=False;
  EdRadPresCoeff.Enabled:=False;
  CoordinateWindow.XcoordInput.Enabled:=True;
  CoordinateWindow.YcoordInput.Enabled:=True;
  CoordinateWindow.ZcoordInput.Enabled:=True;
  CbOrientations.ItemIndex:=0;
  ResultHandler:=MessageDlgEx('Reset','Do you want to redefine the agglomerate or'+
                              ' to alter the parameters such as the Number of '+
                              'Orientations?',mtConfirmation,buttons,results,0,
                              dckActiveForm,0,-1,-1);
  if ResultHandler=1 then
    begin
      for i:=0 to Nl-1 do
        begin
          dia[i]:=0;
          Xcoord[i]:=0;
          Ycoord[i]:=0;
          Zcoord[i]:=0;
          RefrIndexReal[i]:=0;
          RefrIndexImag[i]:=0;
          Rcoord[i]:=0;
          Phicoord[i]:=0;
          Thetacoord[i]:=0;
          RelRefrIndex[i,0]:=0;
          RelRefrIndex[i,1]:=0;
          x[i]:=0;
          WiscombeNumber[i]:=0;
          for j:=0 to Nl-1 do
            begin
              djl[i,j]:=0;
              phijl[i,j]:=0;
              thetajl[i,j]:=0;
            end;
          end;
                //Default settings for global variables
      Nl:=0;
      RefrIndMed:=0;
      WLength:=0;
                //Default settings of the Main Window
      AggrConfigButton.Enabled:=True;
      ResetButton.Enabled:=False;
      CbOrientations.Enabled:=False;
      MainWindow.RefrIndexInput.Enabled:=False;
      MainWindow.WavelengthInput.Enabled:=False;
      MainWindow.StartButton.Enabled:=False;
      MainWindow.ApAngleInput.Text:='';
      MainWindow.ApAngleInput.Enabled:=False;
      MainWindow.Label23.Hide;
    end
  else
    begin
      if length(X0coord)>0 then
        begin
          for i:=0 to Nl-1 do
            begin
              Xcoord[i]:=X0coord[i];
              Ycoord[i]:=Y0coord[i];
              Zcoord[i]:=Z0coord[i];
            end;
        end;
      MainWindow.RefrIndexInput.Enabled:=True;
      MainWindow.RefrIndexInput.SetFocus;
      MainWindow.WavelengthInput.Enabled:=True;
      MainWindow.StartButton.Enabled:=True;
      MainWindow.ApAngleInput.Enabled:=False;
      CbOrientations.Enabled:=True;
      MainWindow.Label23.Show;
      AggrConfigButton.Enabled:=False;
      StartButton.Enabled:=True;
    end;

end;


{-----------------------------------------------------------------------------
  Procedure: TMainWindow.StartButtonClick
  Purpose:   Main Calculation Unit, contains all necessary calculations for
             multisphere scattering
  Comments:  between each process message a delay of 0.5s is included for
             readability
-----------------------------------------------------------------------------}

procedure TMainWindow.StartButtonClick(Sender: TObject);
const
  time=500;
var
  xs,ys,zs:real;
  IL:integer;
  Restarts,i:integer;
  res:real;
  theta_x,theta_y:extended;
  Cext,Csca,AsyPar:array of extended;
  NoO:ShortInt;
  Cext_av,Csca_av,AsyPar_av:extended;
begin
  Label23.Hide;  //"Agglomerate specified!"
  StartButton.Enabled:=False;
  RefrIndexInput.Enabled:=False;
  WaveLengthInput.Enabled:=False;
  CbOrientations.Enabled:=False;
  Label24.Show;  //"Calculation started.."
  Delay(time);
  Label24.Hide;
  Label25.Show;  //"Initializing values.."
  NoO:=StrToInt(CbOrientations.Items[CbOrientations.ItemIndex]);
  SetLength(Cext,NoO);
  SetLength(Csca,NoO);
  SetLength(AsyPar,NoO);
  SetLength(X0coord,length(Xcoord));
  SetLength(Y0coord,length(Ycoord));
  SetLength(Z0coord,length(Zcoord));
  try
        //Reading the values from the Main Window
    RefrIndMed:=StrToFloat(RefrIndexInput.Text);
    WLength:=StrToFloat(WavelengthInput.Text);
    if RefrIndMed<1 then abort;
    if WLength<0 then abort;
        //Wavenumber calculation
    Wavenumber:=1000*TwoPi*RefrIndMed/Wlength;
    Delay(time);
    Label25.Hide;
    Label26.Show;   //"Center of gravity calculation"
        //Geometrical Centre of gravity calculation for the specified agglomerate
    xs:=CenterOfGravity(Xcoord);
    ys:=CenterOfGravity(Ycoord);
    zs:=CenterOfGravity(Zcoord);
    Delay(time);
    Label26.Hide;
    Label27.Show;    //"Point of Origin Transformation"
        //Move origin of coordinate system to Centre of gravity
    PointOfOriginTransformation(xs,ys,zs);
    for i:=0 to Nl-1 do
      begin
        X0coord[i]:=Xcoord[i];
        Y0coord[i]:=Ycoord[i];
        Z0coord[i]:=Zcoord[i];
      end;
    Delay(time);
    Label27.Hide;
    Label33.Show;    //"Transforming to spherical Coords"
        //Calculates Spherical Coordinates of the spheres
    CartesianToSpherical;
    Delay(time);
    Label33.Hide;
    Label28.Show;    //"Calculating Size Parameter"
        //Calculates Size Parameter
    SizeParameter(dia);
    Delay(time);
    Label28.Hide;
    Label29.Show;    //"Calculating Wiscombe criterion"
        //Calculates Wiscombe Criterion for each sphere
    WiscombeCriterion(x);
    Delay(time);
    Label29.Hide;
    Label30.Show;    //"Calculating rel. refractive index"
        //Computes relative refractive index of the spheres
    RelativeRefractiveIndex(RefrIndexReal,RefrIndexImag,RefrIndMed);
    Delay(time);
    Label30.Hide;
    Label32.Show;    //"Calculating Mie Parameter"
        //Computes the Mie Parameters a(n) and b(n)
    Miecalc;
    Delay(time);
    Label32.Hide;
    Label34.Show;   //"Calculating Distance Matrix"
        //Calculates the Distance Matrix of Coordinates between the spheres
    DistanceMatrix;
    Delay(time);
    Label34.Hide;
    Label35.Show;   //"Building Coefficient Matrix"
    Delay(1);
       //Builds the Coefficient Matrix for the calculation of the partial
       //scattering coefficients
    Matrix_Construction;
    Label35.Hide;
    Label36.Show;   //"Solving linear equation system"
    Delay(1);
    SolutionVector:=Copy(Vector_pq,0,length(Vector_pq));
        //Solves the equation system for the partial scattering coefficients
        //a(mn,j), b(mn,j)
    try
      BICGSTAB(Matrix_AB,SolutionVector,1E-6,1E-6,10000,Vector_amn_bmn,IL,
               Restarts,Res);
    except
      on E_ComplexLinearEquationSystemSolver_Error do
        begin
          SolutionVector:=Copy(Vector_pq,0,length(Vector_pq));
          GMRES_m(Matrix_AB,SolutionVector,1E-6,100,1000,Vector_amn_bmn,
                  Restarts,IL,Res);
        end;
    end;
    Label36.Hide;
    SetLength(Matrix_AB,0,0);
    SetLength(SolutionVector,0);
    Label37.Show;   //"Calculating Extinction Coefficient.."
    Delay(time);
    Cext[0]:=ExtinctionCoefficient;
    Label37.Hide;
    if Absorption then
      begin
        Label38.Show;  //"Calculating Scattering Coefficient.."
        Delay(1);
        CalcVector_amnl_bmnl;
        Csca[0]:=ScatteringCoefficient;
        Label38.Hide;
      end
    else
      Csca[0]:=Cext[0];
    Label39.Show;   //"Calculating Asymmetry Parameter.."
    Delay(time);
    if Absorption then
      AsyPar[0]:=AsymmetryParameter(Csca[0])
    else
      begin
        CalcVector_amnl_bmnl;
        AsyPar[0]:=AsymmetryParameter(Csca[0]);
      end;
    Label39.Hide;

        {                      Orientation Averaging                    }

    for i:=1 to NoO-1 do
      begin
        Application.ProcessMessages;
        if CloseClicked then
          exit;
        Label40.Show;   //"Changing Orientation of Agglomerate.."
        Delay(time);
        GetRotationAngles(i,NoO,theta_x,theta_y);
        RotateAgglomerate(X0coord,Y0coord,Z0coord,theta_x,theta_y);
        Label40.Hide;
        Label33.Show;   //"Transforming to Spherical Coords"
        Delay(time);
        CartesianToSpherical;
        Label33.Hide;
        Label34.Show;   //"Calculating Distance Matrix"
        Delay(time);
        DistanceMatrix;
        Label34.Hide;
        Label35.Show;   //"Building Coefficient Matrix"
        Delay(1);
        Matrix_Construction;
        Label35.Hide;
        Label36.Show;   //"Solving linear equation system"
        Delay(1);
        SolutionVector:=Copy(Vector_pq,0,length(Vector_pq));
        GMRES_m(Matrix_AB,SolutionVector,1E-6,100,1000,Vector_amn_bmn,Restarts,IL,Res);
        Label36.Hide;
        SetLength(Matrix_AB,0,0);
        SetLength(SolutionVector,0);
        Label37.Show;   //"Calculating Extinction Coefficient"
        Delay(time);
        Cext[i]:=ExtinctionCoefficient;
        Label37.Hide;
        if Absorption then
          begin
            Label38.Show;       //"Calculating Scattering Coefficient"
            Delay(1);
            CalcVector_amnl_bmnl;
            Csca[i]:=ScatteringCoefficient;
            Label38.Hide;
          end
        else
          Csca[i]:=Cext[i];
        Label39.Show;           //"Calculating Asymmetry Parameter"
        Delay(time);
        if Absorption then
          AsyPar[i]:=AsymmetryParameter(Csca[i])
        else
          begin
            CalcVector_amnl_bmnl;
            AsyPar[i]:=AsymmetryParameter(Csca[i]);
          end;
        Label39.Hide;
      end;

        //Calculating the Average Values of the Cross Sections and the Asymmetry
        //Parameter
    Cext_av:=AverageValue(Cext);
    Csca_av:=AverageValue(Csca);
    AsyPar_av:=AverageValue(AsyPar);
    EdExtCoeff.Text:=FloatToStrF(Cext_av,ffGeneral,6,5);
    EdScatCoeff.Text:=FloatToStrF(Csca_av,ffGeneral,6,5);
    EdAbsCoeff.Text:=FloatToStrF((Cext_av-Csca_av),ffGeneral,6,5);
    EdAsyPar.Text:=FloatToStrF(AsyPar_av,ffGeneral,6,5);
    EdRadPresCoeff.Text:=FloatToStrF(Cext_av-AsyPar_av*Csca_av,ffGeneral,6,5);
    EdExtCoeff.Enabled:=True;
    EdScatCoeff.Enabled:=True;
    EdAbsCoeff.Enabled:=True;
    EdAsyPar.Enabled:=True;
    EdRadPresCoeff.Enabled:=True;

  except
    on EConvertError do
      begin
        ShowMessage(Format(' This is not an allowed value!'+
        '%s All values must be floating points, %s '+
        'Refractive Index must be greater than 1,'+
        ' %s the Wavelength must be positive %s',
        [#10#13,#10#13,#10#13,#10#13,#10#13] ),dckActiveForm);
        RefrIndexInput.SetFocus;
        Label25.Hide;
      end;
  end;
end;


{-----------------------------------------------------------------------------
  Procedure: TMainWindow.FormDestroy
  Purpose:   Realizes the memory release of all global variables if Program
             is shut down
  Comments:
-----------------------------------------------------------------------------}

procedure TMainWindow.FormDestroy(Sender: TObject);
begin
  SetLength(dia,0);
  SetLength(x,0);
  SetLength(Xcoord,0);
  SetLength(Ycoord,0);
  SetLength(Zcoord,0);
  SetLength(Rcoord,0);
  SetLength(Phicoord,0);
  SetLength(Thetacoord,0);
  SetLength(RefrIndexReal,0);
  SetLength(RefrIndexImag,0);
  SetLength(RelRefrIndex,0,0);
  SetLength(MieParameter,0,0,0);
  SetLength(djl,0,0);
  SetLength(phijl,0,0);
  SetLength(thetajl,0,0);
  SetLength(Vector_pq,0);
  SetLength(Matrix_AB,0,0);
  SetLength(Vector_amn_bmn,0);
  SetLength(Vector_amnl_bmnl,0);
  SetLength(WiscombeNumber,0);
end;


{-----------------------------------------------------------------------------
  Procedure: TMainWindow.FormCreate
  Purpose:   Initial Values for the Edit Components of the Main Window
  Comments:
-----------------------------------------------------------------------------}

procedure TMainWindow.FormCreate(Sender: TObject);
begin
  RefrIndexInput.Text:=FloatToStr(1.332);
  WavelengthInput.Text:=FloatToStr(632.8);
end;

end.
