{ plot de paquets d'ondes dans l'espace réel                                    }
{ 16.07.05:FO: première version                                                 }
{ A FAIRE                                                                       }
{ introduire la function d'onde polarisée -> introduire k0+ et k0-              }
{ définir sous la forme d'néergie Zeeman (se rapprocher du SpinEcho)            }
{ calcul 2D pour chaque point de la distribution puis somme des amplitudes      }
{ faire une évaluation du temps de calcul (progress bar et boutton stop         }
{ A FAIRE                                                                       }
{ * c'est l'occasion ou jamais d'utiliser les graphs 3D!!! Youpi                }
{   2D + time                                                                   }
{   l'animation est pas mal non plus                                            }
unit WaveFunctionForm;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, Spin, StdCtrls, EditReal, ExtCtrls, ComCtrls, Menus, Graph2D, ClipBrd,
  math, Math387, GifImage, GR32_Image, ImgList, JvComponent, JvBmpAnim;

const
   NbrMaxPtsDistr = 99;
   rad = 3.14159265358979323846/180;

type
  TWaveFunctionPlotting = class(TForm)
    PageControl1: TPageControl;
    TabSheet1: TTabSheet;
    GroupBox1: TGroupBox;
    Label1: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    Label2: TLabel;
    Theta0Box: TEditReal;
    dThetaBox: TEditReal;
    MultThetaBox: TSpinEdit;
    NbPtsThetaBox: TSpinEdit;
    GroupBox2: TGroupBox;
    Label5: TLabel;
    Label6: TLabel;
    Label7: TLabel;
    Label8: TLabel;
    lambda0box: TEditReal;
    dlambdaBox: TEditReal;
    MultLambdaBox: TSpinEdit;
    NbPtsLambdaBox: TSpinEdit;
    TabSheet2: TTabSheet;
    GroupBox4: TGroupBox;
    Label17: TLabel;
    Label20: TLabel;
    deltaBox: TEditReal;
    Memo1: TMemo;
    Panel2: TPanel;
    status: TLabel;
    ProgressBar: TProgressBar;
    Label18: TLabel;
    ThetaClabel: TLabel;
    TabSheet3: TTabSheet;
    Label9: TLabel;
    XminBox: TEditReal;
    XstepBox: TEditReal;
    Label10: TLabel;
    Label11: TLabel;
    Label15: TLabel;
    X0box: TEditReal;
    XmaxBox: TEditReal;
    Label16: TLabel;
    Label14: TLabel;
    Label13: TLabel;
    Label12: TLabel;
    YminBox: TEditReal;
    YstepBox: TEditReal;
    YmaxBox: TEditReal;
    Y0box: TEditReal;
    MainMenu1: TMainMenu;
    Plot1: TMenuItem;
    STOP1: TMenuItem;
    Graph2D1: TMenuItem;
    Label19: TLabel;
    Label21: TLabel;
    Label22: TLabel;
    Label23: TLabel;
    TabSheet4: TTabSheet;
    Label24: TLabel;
    TminBox: TEditReal;
    Label25: TLabel;
    TstepBox: TEditReal;
    Label26: TLabel;
    TmaxBox: TEditReal;
    Label27: TLabel;
    T0box: TEditReal;
    PlotInVacuumVsT: TMenuItem;
    PlotinvacuuAtT0: TMenuItem;
    Label28: TLabel;
    NeutronSpeedLabel: TLabel;
    Label29: TLabel;
    Save1: TMenuItem;
    SaveDialog: TSaveDialog;
    Animate1: TMenuItem;
    TabSheet5: TTabSheet;
    ImageAnimate: TImage;
    Animation1: TMenuItem;
    ReflectionVsTimenonpolarised1: TMenuItem;
    PacketreflectionondioptreatT01: TMenuItem;
    Label31: TLabel;
    Label30: TLabel;
    Label32: TLabel;
    deltaMbox: TEditReal;
    Label33: TLabel;
    GroupBox5: TGroupBox;
    Label38: TLabel;
    EditReal1: TEditReal;
    Label34: TLabel;
    Label35: TLabel;
    Plotinvacuummagnetic1: TMenuItem;
    Label36: TLabel;
    Label37: TLabel;
    tauNSEbox: TEditReal;
    Label39: TLabel;
    PolarBox: TEditReal;
    Label40: TLabel;
    MangleBox: TEditReal;
    PolarisedNeutronsWithSFbutton: TMenuItem;
    ests1: TMenuItem;
    Vacuum1: TMenuItem;
    Reflectionondioptre1: TMenuItem;
    Magneticreflection1: TMenuItem;
    PhaseBox: TCheckBox;
    SFreflectivityanalytic1: TMenuItem;
    Label41: TLabel;
    GroupBox3: TGroupBox;
    CheckBoxCrop: TCheckBox;
    CheckBoxMerge: TCheckBox;
    CheckBoxPalette: TCheckBox;
    Finitesize1: TMenuItem;
    FiniteSizeTransversecoherenceOnly: TMenuItem;
    Label42: TLabel;
    LengthBox: TEditReal;
    Planewave1: TMenuItem;
    TabSheet6: TTabSheet;
    UnitXbox: TRadioGroup;
    UnitYbox: TRadioGroup;
    GroupBox6: TGroupBox;
    PhiBox: TEditReal;
    dPhiBox: TEditReal;
    Label43: TLabel;
    Label44: TLabel;
    SizeDbox: TEditReal;
    Label45: TLabel;
    dSizeDBox: TEditReal;
    Label46: TLabel;
    OnRandomDomainsButton: TMenuItem;
    Anime1: TMenuItem;
    JvBmpAnimator: TJvBmpAnimator;
    ImageList: TImageList;
    Image1: TImage;
    OnDomainesUPDOWN: TMenuItem;
    GenerateUPDOWNpatternButton: TButton;
    GeneraterandomPatternButton: TButton;
    Label47: TLabel;
    CutPositionXBox: TEditReal;
    Transverse: TMenuItem;
    transvercoherenceonlymagnetic: TMenuItem;
    PLOT1D1: TMenuItem;
    CutInVacuum: TMenuItem;
    CutAfterreflectionondioptre1: TMenuItem;
    CutafterreflectiononmagneticdomainsUPDO1: TMenuItem;
    CutPositionYbox: TEditReal;
    Label48: TLabel;
    Label49: TLabel;
    Label50: TLabel;
    NrDomainsAveragesBox: TSpinEdit;
    Label51: TLabel;
    procedure FormCreate(Sender: TObject);
    procedure deltaBoxKeyDown(Sender: TObject; var Key: Word; Shift: TShiftState);
    procedure STOP1Click(Sender: TObject);
    procedure PlotInVacuumVsTClick(Sender: TObject);
    procedure PlotinvacuuAtT0Click(Sender: TObject);
    procedure lambda0boxKeyDown(Sender: TObject; var Key: Word;Shift: TShiftState);
    procedure Save1Click(Sender: TObject);
    procedure Animate1Click(Sender: TObject);
    procedure ReflectionVsTimenonpolarised1Click(Sender: TObject);
    procedure PacketreflectionondioptreatT01Click(Sender: TObject);
    procedure Plotinvacuummagnetic1Click(Sender: TObject);
    procedure initParameters;
    procedure GenerateDistributions;
    procedure PolarisedNeutronsWithSF(name:string);
    procedure WaveReflection1D(name:string;indice:double);
    procedure WaveReflection(name:string;indice:double);
    procedure WavePacketReflection1D(name:string;indice,T0:double);
    procedure WavePacketReflection(name:string;indice,T0:double);
    procedure SFreflectivityanalytic1Click(Sender: TObject);

    procedure ReflectivityWithSF(k0,thetai,indice,indiceM,phi,T0:double;var ru,rd,rsf:TCplx);
    procedure NeutronSpeedLabelClick(Sender: TObject);
    procedure FiniteSizeTransversecoherenceOnlyClick(Sender: TObject);
    procedure Planewave1Click(Sender: TObject);
    procedure Anime1Click(Sender: TObject);
    procedure Graph2D3DFormStartMenu1Click(Sender: TObject);
    procedure OnDomainesUPDOWNClick(Sender: TObject);
    procedure GenerateUPDOWNpattern(name:string);    
    procedure GenerateUPDOWNpatternButtonClick(Sender: TObject);
    procedure GeneraterandomPattern(name:string);    
    procedure GeneraterandomPatternButtonClick(Sender: TObject);
    procedure TransverseClick(Sender: TObject);
    procedure transvercoherenceonlymagneticClick(Sender: TObject);
    procedure OnDomainsUPDOWN(nameD,nameP:string);
    procedure OnRandomDomainsButtonClick(Sender: TObject);
    procedure CutInVacuumClick(Sender: TObject);
    procedure CutAfterreflectionondioptre1Click(Sender: TObject);
    procedure CutOnDomainsUpDown1D(nameD,nameP:string);
    procedure CutafterreflectiononmagneticdomainsUPDO1Click(Sender: TObject);
    procedure PolarisedNeutronsWithSFbuttonClick(Sender: TObject);
  private
    { Déclarations privées }
  public
    { Déclarations publiques }
  end;

var
  WaveFunctionPlotting: TWaveFunctionPlotting;
  DistrT,DistrTc,DistrTs,DistrK,DistrTint,DistrKint:Array [0..NbrMaxPtsDistr] of double;
  stop:boolean;
  lambda0,dlambda,theta0,dtheta,k0,dk,x0,y0,t0,Xc,Yc:double;
  xmin,xmax,xstep,ymin,ymax,ystep,Xunit,Yunit,tmin,tmax,tstep:double;
  NbPtsK,nbPtsT,MultK,MultT:integer;
  Xsize,Ysize,Tsize,CutSize:integer;
  indice,indiceU,indiceD:double; { indice optique }
  phiM:double;

implementation

uses EditData,MtxVec,constants, Graph2D3D, MainForm, XYcommon;

{$R *.dfm}

procedure TWaveFunctionPlotting.FormCreate(Sender: TObject);
begin
  left:=0;top:=0;width:=335;height:=287;
  WindowState:=wsMinimized;

  include(GIFImageDefaultDrawOptions, goDirectDraw);
end;

procedure TWaveFunctionPlotting.deltaBoxKeyDown(Sender: TObject;var Key: Word; Shift: TShiftState);
var indice,Tc:double;
begin
  { calcul de l'angle d'incidence critique }
  if key=13
  then begin
       indice:=1-deltaBox.value*1e-6;
       Tc:=arcCos(indice)*180/pi;
       ThetaClabel.Caption:=FloatToStrF(Tc,ffFixed,4,6);
       end;
end;

procedure TWaveFunctionPlotting.lambda0boxKeyDown(Sender: TObject;var Key: Word; Shift: TShiftState);
var indice,Tc:double;
begin
  { calcul de la vitesse du neutron }
  if key=13 then begin NeutronSpeedLabel.Caption:=FloatToStrF(hp/mn/(lambda0box.value*1e-9),ffGeneral,3,6); end;
end;

{ calcul de la vitesse pour la longueur d'onde centrale proposée }
procedure TWaveFunctionPlotting.NeutronSpeedLabelClick(Sender: TObject);
begin
  NeutronSpeedLabel.Caption:=FloatToStrF(hp/mn/(lambda0box.value*1e-9),ffGeneral,3,6);
end;


procedure TWaveFunctionPlotting.STOP1Click(Sender: TObject);
begin stop:=true; end;

procedure TWaveFunctionPlotting.initParameters;
begin
   lambda0:=lambda0Box.value;dlambda:=dlambdaBox.value;
   NbPtsK:=NbPtsLambdaBox.value;MultK:=MultLambdaBox.value;
   theta0:=theta0Box.value*rad;dtheta:=dthetaBox.value*rad;
   NbPtsT:=NbPtsThetaBox.value;MultT:=MultThetabox.value;
   k0:=2*Pi/lambda0;dk:=k0*dlambda/lambda0;
   if UnitXbox.itemIndex=0 then Xunit:=1 else Xunit:=1000;
   if UnitYbox.itemIndex=0 then Yunit:=1 else Yunit:=1000;
   xmin:=XminBox.value*Xunit;xstep:=XstepBox.value*Xunit;xmax:=XmaxBox.value*Xunit;Xsize:=round((Xmax-Xmin)/Xstep)+1;
   ymin:=YminBox.value*Yunit;ystep:=YstepBox.value*Yunit;ymax:=YmaxBox.value*Yunit;Ysize:=round((Ymax-Ymin)/Ystep)+1;
   tmin:=TminBox.value;tstep:=TstepBox.value;tmax:=TmaxBox.value;Tsize:=round((Tmax-Tmin)/Tstep)+1;
   x0:=x0box.value*Xunit;y0:=y0box.value*Yunit;t0:=T0box.Value;
   xc:=CutPositionXbox.value*Xunit;yc:=CutPositionYbox.value*Yunit;CutSize:=round(Yc/Ystep);
   indice:=1-deltaBox.value*1e-6; { indice optique }
   indiceU:=1-(deltaBox.value+deltaMbox.value)*1e-6;
   indiceD:=1-(deltaBox.value-deltaMbox.value)*1e-6;
   phiM:=MangleBox.value*pi/180;
end;

{ creation des distributions en theta et en k   }
procedure TWaveFunctionPlotting.GenerateDistributions;
var i,j,nbr:integer;
    sig,delta,totalIntensity:double;
begin
      { generation de la distribution en theta }
      Nbr:=trunc(NbPtsT/2); { 1/2 nbr de points à prendre en compte }
      sig:=dtheta/2; { HWHM }
      totalIntensity:=0;
      for i:=0 to NbPtsT-1 do
        begin
           TotalIntensity:=0;
           delta:=MultT*sig*sqrt(2*ln(2))/NbPtsT; { step dans la distribution }
           for j:=-Nbr to Nbr do { création de la distribution proprement dire }
           begin
              DistrT[j+Nbr]:=theta0+j*delta;
              DistrTc[j+Nbr]:=cos(DistrT[j+Nbr]);
              DistrTs[j+Nbr]:=sin(DistrT[j+Nbr]);
              DistrTInt[j+Nbr]:=1/(sig*sqrt(2*pi))*exp(-(sqr(j*delta)/(2*sig*sig)));
              TotalIntensity:=TotalIntensity+DistrTInt[j+Nbr];  end;
           { renormalisation des intensités } { NECESSAIRE ?! POURQUOI ?! }
           for j:=-Nbr to Nbr do DistrTInt[j+Nbr]:=DistrTInt[j+Nbr]/TotalIntensity;
        end;
      { generation de la distribution en K }
      Nbr:=trunc(NbPtsK/2); { 1/2 nbr de points à prendre en compte }
      sig:=dk/2; { HWHM }
      totalIntensity:=0;
      for i:=0 to NbPtsK-1 do
        begin
           TotalIntensity:=0;
           delta:=MultK*sig*sqrt(2*ln(2))/NbPtsK; { step dans la distribution }
           for j:=-Nbr to Nbr do { création de la distribution proprement dire }
           begin
              DistrK[j+Nbr]:=k0+j*delta;
              DistrKInt[j+Nbr]:=1/(sig*sqrt(2*pi))*exp(-(sqr(j*delta)/(2*sig*sig)));
              TotalIntensity:=TotalIntensity+DistrKInt[j+Nbr];  end;
           { renormalisation des intensités } { NECESSAIRE ?! POURQUOI ?! }
           for j:=-Nbr to Nbr do DistrKInt[j+Nbr]:=DistrKInt[j+Nbr]/TotalIntensity;
        end;
end;

{ ***************************************************************************** }
{ ***************************************************************************** }
{   VACUUM                                                                      }
{ ***************************************************************************** }
{ ***************************************************************************** }

procedure TWaveFunctionPlotting.CutInVacuumClick(Sender: TObject);
var i,k,m,n:integer;
    NbT,NbK:integer;
    indMod,indArg:integer;
    auxX,auxY,auxR1,auxR2:double;
    x,y,phase,w,T,wT:double;
begin
   InitParameters;generateDistributions;
   indMod:=und.addspectre('Phi_mod',CutSize,1,1);
   indArg:=und.addspectre('Phi_arg',CutSize,1,1);
   NbT:=trunc(NbPtsT/2);NbK:=trunc(NbPtsK/2);
   for i:=0 to CutSize-1 do
   begin
      X:=Xc;Y:=i*Ystep;
      auxX:=0;auxY:=0;
      for m:=-NbT to NbT do
        for n:=-NbK to NbK do
        begin
          auxR1:=DistrTint[m+NbT]*DistrKint[n+NbK];
          auxR2:=DistrK[n+NbK]*(DistrTc[m+NbT]*(X-X0)+DistrTs[m+NbT]*(Y-Y0));
          auxX:=auxX+AuxR1*cos(auxR2-wt);
          auxY:=auxY-AuxR1*sin(auxR2-wt);
        end;
      und.Spect[indMod].data[i,0,0].x:=Y/Yunit;
      und.Spect[indMod].data[i,0,0].z:=sqr(auxX)+sqr(auxY); {Cabs(aux); }
      und.Spect[indArg].data[i,0,0].x:=Y/Yunit;
      und.Spect[indArg].data[i,0,0].z:=Arg(Cplx(auxX,auxY))/rad;
   end;
end;     { fin Plotting Vs T0 }

{ plot de la fonction d'onde à T=T0                                             }
{ toutes les unités sont en nm                                                  }
{ on va utiliser des distributions gaussiennes                                  }
{ semble fonctionner normalement                                                }
procedure TWaveFunctionPlotting.PlotinvacuuAtT0Click(Sender: TObject);
var i,j,k,m,n:integer;
    NbT,NbK:integer;
    indMod,indArg:integer;
    auxX,auxY,auxR1,auxR2:double;
    x,y,phase,w,T,wT:double;
begin
   stop:=false;
   status.Caption:='Calculating...';status.Update;
   ProgressBar.max:=Xsize;ProgressBar.position:=0;ProgressBar.update;

   InitParameters;generateDistributions;
   indMod:=deuxd.addspectre('Phi_mod',Xsize,Ysize,1);
   indArg:=deuxd.addspectre('Phi_arg',Xsize,Ysize,1);

   { on calcul la phase sans dispersion en k, sinon le paquet d'onde est trop étalé }
   wt:=hb*sqr(k0*1e9)*(T0*1e-9)/mn; { E = hb w / 2 ; on fait le calcul en SI }
   { methode pedestre ; point par point dans le plan (X,Y) }
   { plus efficace que la méthode vectorielle !!!          }
   NbT:=trunc(NbPtsT/2);
   NbK:=trunc(NbPtsK/2);
   for i:=0 to Xsize-1 do
   begin
     ProgressBar.position:=i;ProgressBar.update;
     for j:=0 to Ysize-1 do
     begin
       X:=xmin+i*Xstep;Y:=ymin+j*Ystep;
       auxX:=0;auxY:=0;
       for m:=-NbT to NbT do
         for n:=-NbK to NbK do
         begin
{            wt:=hb*sqr(DistrK[n+NbK]*1e9)*(T0*1e-9)/mn/2; { on fait le calcul en SI }
            auxR1:=DistrTint[m+NbT]*DistrKint[n+NbK];
            auxR2:=DistrK[n+NbK]*(DistrTc[m+NbT]*(X-X0)+DistrTs[m+NbT]*(Y-Y0));
            auxX:=auxX+AuxR1*cos(auxR2-wt);
            auxY:=auxY-AuxR1*sin(auxR2-wt);
         end;
       deuxd.Spect[indMod].data[i,j,0].x:=X/Xunit;
       deuxd.Spect[indMod].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indMod].data[i,j,0].z:=sqr(auxX)+sqr(auxY); {Cabs(aux); }
       deuxd.Spect[indArg].data[i,j,0].x:=X/Xunit;
       deuxd.Spect[indArg].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indArg].data[i,j,0].z:=Arg(Cplx(auxX,auxY))/rad;
     end;
     Application.processMessages;if Stop then break;
   end;
   status.Caption:='Idle';status.Update;
end;     { fin Plotting Vs T0 }


{ plot de la fonction d'onde dans l'espace réel en fonction du temps            }
{ toutes les unités sont en nm                                                  }
{ on va utiliser des distributions gaussiennes                                  }
procedure TWaveFunctionPlotting.PlotInVacuumVsTClick(Sender: TObject);
var i,j,k,m,n:integer;
    NbT,NbK:integer;
    indMod,indArg:integer;
    X,Y,auxX,auxY,auxR1,auxR2:double;
    phase,w,T,wT:double;
begin
   stop:=false;
   status.Caption:='Calculating...';status.Update;
   initParameters;generateDistributions;
   ProgressBar.max:=Tsize;ProgressBar.position:=0;

   for k:=1 to Tsize do
   begin
   indMod:=deuxd.addspectre('Phi_mod'+IntToStr(k),Xsize,Ysize,1);
   { indArg:=deuxd.addspectre('Phi_arg',Xsize,Ysize,1); }
   { calcul de omega } { w = hbar k² / m }
   T:=Tmin+(k-1)*Tstep;
   NbT:=trunc(NbPtsT/2);NbK:=trunc(NbPtsK/2);
   for i:=0 to Xsize-1 do
   begin
     for j:=0 to Ysize-1 do
     begin
       X:=xmin+i*Xstep;Y:=ymin+j*Ystep;
       auxX:=0;auxY:=0;
       for m:=-NbT to NbT do
         for n:=-NbK to NbK do
         begin
            wt:=hb*sqr(DistrK[n+NbK])*T/mn/2;
            auxR1:=DistrTint[m+NbT]*DistrKint[n+NbK];
            auxR2:=DistrK[n+NbK]*(DistrTc[m+NbT]*(X-X0)+DistrTs[m+NbT]*(Y-Y0));
            auxX:=auxX+AuxR1*cos(auxR2-wt);
            auxY:=auxY-AuxR1*sin(auxR2-wt);
         end;
       deuxd.Spect[indMod].data[i,j,0].x:=X/Xunit;
       deuxd.Spect[indMod].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indMod].data[i,j,0].z:=sqr(auxX)+sqr(auxY); {Cabs(aux); }
{       deuxd.Spect[indArg].data[i,j,0].x:=X;
       deuxd.Spect[indArg].data[i,j,0].y:=Y;
       deuxd.Spect[indArg].data[i,j,0].z:=Arg(Cplx(auxX,auxY))/rad; }
     end;
     Application.processMessages;if Stop then break;
     ProgressBar.position:=k;ProgressBar.update;
   end;
   end;
   status.Caption:='Idle';status.Update;
end;     { fin Plotting Vs Time }

{ ***************************************************************************** }
{ ***************************************************************************** }
{     DIOPTRE                                                                   }
{ ***************************************************************************** }
{ ***************************************************************************** }

{ calcul de la fonction d'onde au point (X,Y)                                   }
{ INPUT : la position dans l'espace direct (X,Y)                                }
{         le vecteur d'onde et la dircetion de l'onde plane                     }
{ OUTPUT: les composantes réelles imaginaire de l'onde plane incidente et réfléchie }
procedure calculPsiR(X,Y,p,q,k1,phase:double;var PsiX,PsiY:double);
var q1,R,Rx,Ry,aux,auxI,auxR:double;
begin
        if (k1>=p)
        then begin { réflexion partielle } { ok }
             q1:=sqrt(sqr(k1)-sqr(p)); { vecteur d'onde réel }
             R:=(q-q1)/(q+q1);  { coefficient de réflexion réel ; pas de facteur de phase...curieux... }
             auxI:=p*(X-X0)-q*(Y-Y0)-phase;             { facteur de phase de l'onde incidente }
             auxR:=p*(X-X0)+q*(Y-Y0)-phase;             { facteur de phase de l'onde réfléchie }
             psiX:=cos(auxI)+R*cos(auxR);
             psiY:=sin(auxI)+R*sin(auxR); end
        else begin { onde totalement réfléchie } {ok }
             q1:=sqrt(-sqr(k1)+sqr(p)); { composante du vecteur d'onde imaginaire pure }
             aux:=(sqr(q)+sqr(q1));
             Rx:=(sqr(q)-sqr(q1))/aux;Ry:=-2*q*q1/aux;
             auxI:=p*(X-X0)-q*(Y-Y0)-phase;             { facteur de phase de l'onde incidente }
             auxR:=p*(X-X0)+q*(Y-Y0)-phase;
             psiX:=cos(auxI)+Rx*cos(auxR)-Ry*sin(auxR);
             psiY:=sin(auxI)+Rx*sin(auxR)+Ry*cos(auxR); end
end;

{ calcul de la fonction d'onde au point (X,Y)                                   }
{ INPUT : la position dans l'espace direct (X,Y)                                }
{         le vecteur d'onde et la dircetion de l'onde plane                     }
{ OUTPUT: les composantes réelles imaginaire de l'onde plane incidente et réfléchie }
procedure calculPsiT(X,Y,p,q,k1,phase:double;var PsiX,PsiY:double);
var q1,T,aux,auxT:double;
begin
   if (k1>=p)
   then begin { réflexion partielle } { ok }
        q1:=sqrt(+sqr(k1)-sqr(p));
        T:=2*q/(q+q1); { coefficient de transmission }
        auxT:=p*(X-X0)-q1*(Y-Y0)-phase;
        PsiX:=T*cos(auxT);
        PsiY:=T*sin(auxT); end
   else begin { onde evanescente } { ok }
        q1:=sqrt(-sqr(k1)+sqr(p)); { composantes imaginaire pure ! }
        aux:=2*q/(sqr(q)+sqr(q1))*exp(q1*y);
        auxT:=p*(X-X0)-phase;
        PsiX:=aux*(q*cos(auxT)+q1*sin(auxT));
        PsiY:=aux*(q*sin(auxT)-q1*cos(auxT)); end;
end;

{ plot de la fonction d'onde se reflechissant sur un dioptre                    }
{ INPUT : name, name of the spectrum into which the result must be stored       }
{         indice, optical index                                                 }
procedure TWaveFunctionPlotting.WaveReflection1D(name:string;indice:double);
var i,m,n,NbT:integer;
    indMod,indArg:integer;
    X,Y,psiX,PsiY,auxX,auxY,intT:double;
    p,q,k1:double;
    PhaseYN:boolean;
begin
   indMod:=und.addspectre(name+'_mod',Cutsize,1,1);
   PhaseYN:=PhaseBox.checked;if PhaseYN then indArg:=und.addspectre(name+'_arg',Cutsize,1,1);
   NbT:=trunc(NbPtsT/2);
   for i:=0 to Cutsize-1 do
   begin
       X:=Xc;Y:=i*Ystep;
       auxX:=0;auxY:=0;
       for m:=-NbT to NbT do
       begin
         intT:=DistrTint[m+NbT];
         { calcul de la fonction d'onde }
         p:=k0*DistrTc[m+NbT]; { composante k //    }
         q:=k0*DistrTs[m+NbT]; { composante k perp  }
         k1:=indice*k0;        { vecteur d'onde dans le dioptre     }
         if (Y>0)
         then begin { ********** au-dessus du dioptre  ********************  }
              calculPsiR(X,Y,p,q,k1,0,PsiX,PsiY);
              auxX:=auxX+intT*PsiX;auxY:=auxY+intT*PsiY; end
         else begin { ********** en-dessous du dioptre  ********************  }
              calculPsiT(X,Y,p,q,k1,0,PsiX,PsiY);
              auxX:=auxX+intT*PsiX;auxY:=auxY+intT*PsiY; end
       end;
       und.Spect[indMod].data[i,0,0].x:=Y/Yunit;
       und.Spect[indMod].data[i,0,0].z:=sqr(auxX)+sqr(auxY);
       if PhaseYN then begin
       und.Spect[indArg].data[i,0,0].x:=Y/Yunit;
       und.Spect[indArg].data[i,0,0].z:=Arg(Cplx(auxX,auxY))/rad; end;
   end;
end;

procedure TWaveFunctionPlotting.Cutafterreflectionondioptre1Click(Sender: TObject);
begin
   initParameters;generateDistributions;
   WaveReflection1D('Phi',indice);
   sendmessage(principalform.handle,wm_redessine,0,0);
end;

{ plot de la fonction d'onde se reflechissant sur un dioptre                    }
{ INPUT : name, name of the spectrum into which the result must be stored       }
{         indice, optical index                                                 }
procedure TWaveFunctionPlotting.WaveReflection(name:string;indice:double);
var i,j,m,n,NbT:integer;
    indMod,indArg:integer;
    X,Y,psiX,PsiY,auxX,auxY,intT:double;
    p,q,k1:double;
    PhaseYN:boolean;
begin
   stop:=false;
   status.Caption:='Calculating...';status.Update;
   indMod:=deuxd.addspectre(name+'_mod',Xsize,Ysize,1);
   PhaseYN:=PhaseBox.checked;if PhaseYN then indArg:=deuxd.addspectre(name+'_arg',Xsize,Ysize,1);
   ProgressBar.max:=Xsize-1;ProgressBar.position:=0;progressBar.update;
   NbT:=trunc(NbPtsT/2);
   for i:=0 to Xsize-1 do
   begin
     for j:=0 to Ysize-1 do
     begin
       X:=xmin+i*Xstep;Y:=ymin+j*Ystep;
       auxX:=0;auxY:=0;
       for m:=-NbT to NbT do
       begin
         intT:=DistrTint[m+NbT];
         { calcul de la fonction d'onde }
         p:=k0*DistrTc[m+NbT]; { composante k //    }
         q:=k0*DistrTs[m+NbT]; { composante k perp  }
         k1:=indice*k0;        { vecteur d'onde dans le dioptre     }
         if (Y>0)
         then begin { ********** au-dessus du dioptre  ********************  }
              calculPsiR(X,Y,p,q,k1,0,PsiX,PsiY);
              auxX:=auxX+intT*PsiX;auxY:=auxY+intT*PsiY; end
         else begin { ********** en-dessous du dioptre  ********************  }
              calculPsiT(X,Y,p,q,k1,0,PsiX,PsiY);
              auxX:=auxX+intT*PsiX;auxY:=auxY+intT*PsiY; end
       end;
       deuxd.Spect[indMod].data[i,j,0].x:=X/Xunit;deuxd.Spect[indMod].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indMod].data[i,j,0].z:=sqr(auxX)+sqr(auxY);
       if PhaseYN then begin
       deuxd.Spect[indArg].data[i,j,0].x:=X/Xunit;deuxd.Spect[indArg].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indArg].data[i,j,0].z:=Arg(Cplx(auxX,auxY))/rad; end;
     end;
     Application.processMessages;if Stop then break;
     ProgressBar.position:=i;ProgressBar.update;
   end;
   status.Caption:='Idle';status.Update;
end;


{ plot de la fonction d'onde se reflechissant sur un dioptre                    }
{ à l'instant T0                                                                }
{ INPUT : name, name of the spectrum into which the result must be stored       }
{         indice, optical index                                                 }
{ ATTENTION : indice redevient une variable locale !!!                          }
procedure TWaveFunctionPlotting.WavePacketReflection1D(name:string;indice,T0:double);
var i,m,n:integer;
    NbT,NbK:integer;
    indMod,indArg:integer;
    X,Y,psiX,PsiY,auxX,auxY:double;
    p,q,k1:double;
    intK,intT,wt:double;
    PhaseYN:boolean;
begin
   indMod:=und.addspectre(name+'_mod',CutSize,1,1);
   PhaseYN:=PhaseBox.checked;if PhaseYN then indArg:=und.addspectre(name+'_arg',CutSize,1,1);
   NbT:=trunc(NbPtsT/2);NbK:=trunc(NbPtsK/2);
   for i:=0 to CutSize-1 do
   begin
      X:=Xc;Y:=i*Ystep;
      auxX:=0;auxY:=0;
      for m:=-NbT to NbT do
        for n:=-NbK to NbK do
         begin
            intK:=DistrKint[n+NbK];
            intT:=DistrTint[m+NbT];
            k0:=DistrK[n+NbK];
            wt:=hb*sqr(k0)*T0/mn/2;
            { calcul de la fonction d'onde }
            p:=k0*DistrTc[m+NbT]; { composante k //    }
            q:=k0*DistrTs[m+NbT]; { composante k perp  }
            k1:=indice*k0;        { vecteur d'onde dans le dioptre     }
            if (Y>0)
            then begin { ********** au-dessus du dioptre  ********************  }
                 calculPsiR(X,Y,p,q,k1,wt,PsiX,PsiY);
                 auxX:=auxX+intK*intT*PsiX;auxY:=auxY+intK*intT*PsiY; end
            else begin { ********** en-dessous du dioptre  ********************  }
                 calculPsiT(X,Y,p,q,k1,wt,PsiX,PsiY);
                 auxX:=auxX+intK*intT*PsiX;auxY:=auxY+intK*intT*PsiY; end
         end;
       und.Spect[indMod].data[i,0,0].x:=X/Xunit;
       und.Spect[indMod].data[i,0,0].z:=sqr(auxX)+sqr(auxY);
       if PhaseYN then begin
       und.Spect[indArg].data[i,0,0].x:=X/Xunit;
       und.Spect[indArg].data[i,0,0].z:=Arg(Cplx(auxX,auxY))/rad; end;
   end;
end;

{ plot de la fonction d'onde se reflechissant sur un dioptre                    }
{ à l'instant T0                                                                }
{ INPUT : name, name of the spectrum into which the result must be stored       }
{         indice, optical index                                                 }
{ ATTENTION : indice redevient une variable locale !!!                          }
procedure TWaveFunctionPlotting.WavePacketReflection(name:string;indice,T0:double);
var i,j,m,n:integer;
    NbT,NbK:integer;
    indMod,indArg:integer;
    X,Y,psiX,PsiY,auxX,auxY:double;
    p,q,k1:double;
    intK,intT,wt:double;
    PhaseYN:boolean;
begin
   stop:=false;
   status.Caption:='Calculating...';status.Update;
   indMod:=deuxd.addspectre(name+'_mod',Xsize,Ysize,1);
   PhaseYN:=PhaseBox.checked;if PhaseYN then indArg:=deuxd.addspectre(name+'_arg',Xsize,Ysize,1);
   ProgressBar.max:=Xsize-1;ProgressBar.position:=0;progressBar.update;
   NbT:=trunc(NbPtsT/2);
   NbK:=trunc(NbPtsK/2);
   for i:=0 to Xsize-1 do
   begin
     for j:=0 to Ysize-1 do
     begin
       X:=xmin+i*Xstep;Y:=ymin+j*Ystep;
       auxX:=0;auxY:=0;
       for m:=-NbT to NbT do
         for n:=-NbK to NbK do
         begin
            intK:=DistrKint[n+NbK];
            intT:=DistrTint[m+NbT];
            k0:=DistrK[n+NbK];
            wt:=hb*sqr(k0)*T0/mn/2;
            { calcul de la fonction d'onde }
            p:=k0*DistrTc[m+NbT]; { composante k //    }
            q:=k0*DistrTs[m+NbT]; { composante k perp  }
            k1:=indice*k0;        { vecteur d'onde dans le dioptre     }
            if (Y>0)
            then begin { ********** au-dessus du dioptre  ********************  }
                 calculPsiR(X,Y,p,q,k1,wt,PsiX,PsiY);
                 auxX:=auxX+intK*intT*PsiX;auxY:=auxY+intK*intT*PsiY; end
            else begin { ********** en-dessous du dioptre  ********************  }
                 calculPsiT(X,Y,p,q,k1,wt,PsiX,PsiY);
                 auxX:=auxX+intK*intT*PsiX;auxY:=auxY+intK*intT*PsiY; end
         end;
       deuxd.Spect[indMod].data[i,j,0].x:=X/Xunit;deuxd.Spect[indMod].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indMod].data[i,j,0].z:=sqr(auxX)+sqr(auxY);
       if PhaseYN then begin
       deuxd.Spect[indArg].data[i,j,0].x:=X/Xunit;deuxd.Spect[indArg].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indArg].data[i,j,0].z:=Arg(Cplx(auxX,auxY))/rad; end;
     end;
     Application.processMessages;if Stop then break;
     ProgressBar.position:=i;ProgressBar.update;
   end;
   status.Caption:='Idle';status.Update;
end;

procedure TWaveFunctionPlotting.TransverseClick(Sender: TObject);
begin
   initParameters;generateDistributions;
   WaveReflection('Phi',indice);
   sendmessage(principalform.handle,wm_redessine,0,0);
end;

{ plot de la fonction d'onde se reflechissant sur un dioptre                    }
{ à l'instant T0                                                                }
procedure TWaveFunctionPlotting.PacketreflectionondioptreatT01Click(Sender: TObject);
begin
   initParameters;generateDistributions;
   WavePacketReflection('Phi',indice,T0);
   sendmessage(principalform.handle,wm_redessine,0,0);
end;

{ plot de la fonction d'onde se reflechissant sur un dioptre en fonction du temps }
procedure TWaveFunctionPlotting.ReflectionVsTimenonpolarised1Click(Sender: TObject);
var k:integer;
    instant:double;
begin
  initParameters;generateDistributions;
  for k:=1 to Tsize do
  begin
    instant:=Tmin+(k-1)*Tstep;
    WavePacketReflection('Phi'+IntToStr(k),indice,instant);
  end;
end;

{ ***************************************************************************** }
{ ***************************************************************************** }
{   MAGNETIC CASE                                                               }
{ ***************************************************************************** }
{ ***************************************************************************** }

procedure TWaveFunctionPlotting.ReflectivityWithSF(k0,thetai,indice,indiceM,phi,T0:double;var ru,rd,rsf:TCplx);
var i,indUU,indDD,indSF,size:integer;
    indiceU,indiceD:double;
    k1,k1u,k1d,Q,p0,Ti:double;
    cphi,sphi:double;
    q0,q1D,q1U,c2u,s2u,c2d,s2d,su,sd,du,dd:TCplx;
begin
   indiceU:=indice+indiceM;indiceD:=indice-indiceM;
   cphi:=cos(phi/2);sphi:=sin(phi/2);
   c2u:=Cplx(sqr(cphi),0);s2u:=Cplx(sqr(sphi),0);
   c2d:=Cplx(sqr(sphi),0);s2d:=Cplx(sqr(cphi),0);
   k1:=k0*indice;k1U:=k0*indiceU;k1D:=k0*indiceD;
   p0:=k0*cos(Ti);  { projection du vecteur d'onde sur (Ox) }
   q0:=Cplx(k0*sin(Ti),0);
   q1u:=Csqrt(Cplx(sqr(k1U)-sqr(p0),0));
   q1d:=Csqrt(Cplx(sqr(k1D)-sqr(p0),0));
   su:=Cadd(q0,q1u);sd:=Cadd(q0,q1d); { somme UP ; somme DOWN }
   du:=Csub(q0,q1u);dd:=Csub(q0,q1d); { diff UP  ; diff DOWN }
   ru:=Cdiv(Cadd( Cmul(c2u,Cmul(du,sd)),Cmul(s2u,Cmul(dd,su))),Cadd( Cmul(c2u,Cmul(su,sd)),Cmul(s2u,Cmul(sd,su))) );
   rd:=Cdiv(Cadd( Cmul(c2d,Cmul(du,sd)),Cmul(s2d,Cmul(dd,su))),Cadd( Cmul(c2d,Cmul(su,sd)),Cmul(s2d,Cmul(sd,su))) );
   rsf:=Cmul(Cmul(Cplx(2*cphi*sphi,0),q0),Cdiv( Csub(q1u,q1d), Cadd( Cmul(c2u,Cmul(su,sd)),Cmul(s2u,Cmul(sd,su)))));
end;

{ calcul de la reflectivité SF avec des formules analytiques                    }
procedure TWaveFunctionPlotting.SFreflectivityanalytic1Click(Sender: TObject);
const qmax=0.5;
      qstep=0.002;  { en nm-1 }
      lambda=1;         { lambda = 1nm par défaut }
var i,indUU,indDD,indSF,size:integer;
    indice,indiceU,indiceD,phi:double;
    k1,k1u,k1d,Q,p0,Ti,Ru,Rd,Rsf:double;
    cphi,sphi:double;
    q0,q1D,q1U,c2u,s2u,c2d,s2d,su,sd,du,dd:TCplx;
begin
   indice:=1-(deltaBox.value)*1e-6;
   indiceU:=1-(deltaBox.value+deltaMbox.value)*1e-6;
   indiceD:=1-(deltaBox.value-deltaMbox.value)*1e-6;
   phi:=MangleBox.value*pi/180;
   cphi:=cos(phi/2);sphi:=sin(phi/2);
   c2u:=Cplx(sqr(cphi),0);s2u:=Cplx(sqr(sphi),0);
   c2d:=Cplx(sqr(sphi),0);s2d:=Cplx(sqr(cphi),0);
   size:=round(Qmax/Qstep);
   indUU:=und.addspectre('Ruu',size,1,1);
   indDD:=und.addspectre('Rdd',size,1,1);
   indSF:=und.addspectre('Rsf',size,1,1);
   k0:=2*pi/lambda;
   k1:=k0*indice;k1U:=k0*indiceU;k1D:=k0*indiceD;
   for i:=0 to size-1 do
   begin
      Q:=(i+1)*qstep; { composante du vecteur d'onde perpendiculaire }
      und.spect[indUU].data[i,0,0].x:=Q;
      und.spect[indDD].data[i,0,0].x:=Q;
      und.spect[indSF].data[i,0,0].x:=Q;
      Ti:=arcsin(Q*lambda/4/pi);
      p0:=k0*cos(Ti);  { projection du vecteur d'onde sur (Ox) }
      q0:=Cplx(k0*sin(Ti),0);
      q1u:=Csqrt(Cplx(sqr(k1U)-sqr(p0),0));
      q1d:=Csqrt(Cplx(sqr(k1D)-sqr(p0),0));
      su:=Cadd(q0,q1u);sd:=Cadd(q0,q1d); { somme UP ; somme DOWN }
      du:=Csub(q0,q1u);dd:=Csub(q0,q1d); { diff UP  ; diff DOWN }
      Ru:=sqr(Cabs(Cdiv(Cadd( Cmul(c2u,Cmul(du,sd)),
                              Cmul(s2u,Cmul(dd,su))),
                        Cadd( Cmul(c2u,Cmul(su,sd)),
                              Cmul(s2u,Cmul(sd,su))) )));
      Rd:=sqr(Cabs(Cdiv(Cadd( Cmul(c2d,Cmul(du,sd)),
                              Cmul(s2d,Cmul(dd,su))),
                        Cadd( Cmul(c2d,Cmul(su,sd)),
                              Cmul(s2d,Cmul(sd,su))) )));
      Rsf:=sqr(2*Cabs(q0)*cphi*sphi*
               Cabs( Cdiv( Csub(q1u,q1d), Cadd( Cmul(c2u,Cmul(su,sd)),
                                                Cmul(s2u,Cmul(sd,su))) )));

      und.spect[indUU].data[i,0,0].z:=Ru;
      und.spect[indDD].data[i,0,0].z:=Rd;
      und.spect[indSF].data[i,0,0].z:=Rsf;
   end;
end;

{ calcul des fonctions d'ondes pour polar UP et polar DOWN                      }
procedure TWaveFunctionPlotting.transvercoherenceonlymagneticClick(Sender: TObject);
begin
   initParameters;generateDistributions;
   WaveReflection('PhiUP',indiceU);
   WaveReflection('PhiDO',indiceD);
end;

{ plot de la fonction d'onde se reflechissant sur un dioptre en fonction du temps }
{ calcul séparé pour UP et DOWN -> on considère smplement les indices UP et DOWN  }
{ on introduit un offset tNSE entre UP et DOWN                                    }
{ OUTPUT : les modules pour UP et DOWN (2 spectres incohérents)                   }
{ A FAIRE: somme cohérente des deux ondes                                         }
{ il faudrait rajouter k+ et k-                                                   }
{ IMO la différence est négligeable au niveau de l'échantillon                    }
procedure TWaveFunctionPlotting.Plotinvacuummagnetic1Click(Sender: TObject);
var tauNSE:double;
begin
   initParameters;generateDistributions;
   tauNSE:=tauNSEbox.value; { SE time in a.u., same at time }
   WavePacketReflection('PhiUP',indiceU,T0+tauNSE);
   WavePacketReflection('PhiDO',indiceD,T0-tauNSE);
end;

{ calcul de la fonction d'onde au point (X,Y)                                   }
{ INPUT : la position dans l'espace direct (X,Y)                                }
{         le vecteur d'onde et la dircetion de l'onde plane                     }
{ OUTPUT: les composantes réelles imaginaire de l'onde plane incidente et réfléchie }
procedure calculPsiRpol(X,Y,p0,q,k1U,k1D,phiM,Psi0Ur,Psi0Ui,Psi0Dr,Psi0Di:double;
                        var PsiUr,PsiUi,PsiDr,PsiDi:double);
var cphi,sphi:double;
    q0,q1u,q1d,su,du,sd,dd,c2u,s2u,c2d,s2d:TCplx;
    RU,RD,RSF:TCplx;
    PhaseIr,PhaseIi,PhaseRr,PhaseRi,PsiRUr,PsiRUi,PsiRDr,PsiRDi:double;
begin
   cphi:=cos(phiM/2);sphi:=sin(phiM/2);
   c2u:=Cplx(sqr(cphi),0);s2u:=Cplx(sqr(sphi),0);
   c2d:=Cplx(sqr(sphi),0);s2d:=Cplx(sqr(cphi),0);
   q0:=Cplx(q,0);
   q1u:=Csqrt(Cplx(sqr(k1U)-sqr(p0),0));
   q1d:=Csqrt(Cplx(sqr(k1D)-sqr(p0),0));
   su:=Cadd(q0,q1u);sd:=Cadd(q0,q1d); { somme UP ; somme DOWN }
   du:=Csub(q0,q1u);dd:=Csub(q0,q1d); { diff UP  ; diff DOWN }
   RU:=Cdiv(Cadd( Cmul(c2u,Cmul(du,sd)), Cmul(s2u,Cmul(dd,su))),
            Cadd( Cmul(c2u,Cmul(su,sd)), Cmul(s2u,Cmul(sd,su))) ); { ok }
   RD:=Cdiv(Cadd( Cmul(c2d,Cmul(du,sd)), Cmul(s2d,Cmul(dd,su))),
            Cadd( Cmul(c2d,Cmul(su,sd)), Cmul(s2d,Cmul(sd,su))) ); { ok }
   RSF:=Cmul(Cplx(2*q*cphi*sphi,0),
        Cdiv( Csub(q1u,q1d), Cadd( Cmul(c2u,Cmul(su,sd)),
                                   Cmul(s2u,Cmul(sd,su))) ));  
   { on ne fait pas la différence entre onde totalement réfléchie et onde       }
   { partiellement réfléchie, cela n'a plus de sens                             }
   PhaseIr:=cos(p0*(X-X0)-q*(Y-Y0));PhaseIi:=sin(p0*(X-X0)-q*(Y-Y0)); { facteur de phase de l'onde incidente }
   PhaseRr:=cos(p0*(X-X0)+q*(Y-Y0));PhaseRi:=sin(p0*(X-X0)+q*(Y-Y0)); { facteur de phase de l'onde réfléchie }

   psiRUr:= Psi0Ur*RU.re  + Psi0Dr*RSF.re;
   psiRUi:= Psi0Ur*RU.im  + Psi0Dr*RSF.im;
   psiRDr:= Psi0Ur*RSF.re + Psi0Dr*RD.re;
   psiRDi:= Psi0Ur*RSF.im + Psi0Dr*RD.im;

   psiUr:= psi0Ur*phaseIr + psiRUr*phaseRr - psiRUi*phaseRi;
   psiUi:= psi0Ur*phaseIi + psiRUr*phaseRi + psiRUi*phaseRr;
   psiDr:= psi0Dr*phaseIr + psiRDr*phaseRr - psiRDi*phaseRi;
   psiDi:= psi0Dr*phaseIi + psiRDr*phaseRi + psiRDi*phaseRr;
end;

{ calcul de la reflectivité sur une couche magnétique d'orientation quelconque  }
{ OUTPUT: |Psi+|² et |Psi-|²                                                    }
{         la polarisation Pp Pm                                                 }
{ Pp et Pm désignent les composantes |+> et |->                                 }
{ plot de la polariation avec un code couleur : rouge +, blanc, non pol, bleu - }
procedure TWaveFunctionPlotting.PolarisedNeutronsWithSF(name:string);
var polar:double;
    i,j,m,n:integer;
    NbT,NbK:integer;
    indModU,indModD,indArgU:integer;
    X,Y,psiX,PsiY,auxUX,auxUY,auxDX,auxDY:double;
    p,q,k1p,k1m:double;
    intK,intT,wt,Pp,Pm:double;
    PhaseYN:boolean;
    PsiUr,PsiUi,PsiDr,PsiDi:double;
begin
  { composantes du spineur incident }
  polar:=polarBox.value;
  if polar=1
  then begin Pp:=1;Pm:=0 end
  else begin Pm:=1;Pp:=(1+polar)/(1-polar)*Pm end;
  { pour le moment on ne normalisa pas l'amplitude de l'onde incidente }
  stop:=false;
  status.Caption:='Calculating...';status.Update;
  indModU:=deuxd.addspectre(name+'UP_mod',Xsize,Ysize,1);
  indModD:=deuxd.addspectre(name+'DO_mod',Xsize,Ysize,1);
  PhaseYN:=PhaseBox.checked;if PhaseYN then indArgU:=deuxd.addspectre(name+'U_arg',Xsize,Ysize,1);
  ProgressBar.max:=Xsize-1;ProgressBar.position:=0;progressBar.update;
  NbT:=trunc(NbPtsT/2);
  for i:=0 to Xsize-1 do
  begin
    for j:=0 to Ysize-1 do
    begin
      X:=xmin+i*Xstep;Y:=ymin+j*Ystep;
      auxUX:=0;auxUY:=0;auxDX:=0;auxDY:=0;
      for m:=-NbT to NbT do
      begin
         intT:=DistrTint[m+NbT];
         { calcul de la fonction d'onde }
         p:=k0*DistrTc[m+NbT]; { composante k //    }
         q:=k0*DistrTs[m+NbT]; { composante k perp  }
         k1p:=indiceU*k0;k1m:=indiceD*k0;  { vecteurs d'onde dans le dioptre     }
         if (Y>0)
         then begin { ********** au-dessus du dioptre  ********************  }
              calculPsiRpol(X,Y,p,q,k1p,k1m,phiM,Pp,0,Pm,0,
                            PsiUr,PsiUi,PsiDr,PsiDi);
              auxUX:=auxUX+intT*PsiUr;auxUY:=auxUY+intT*PsiUi;
              auxDX:=auxDX+intT*PsiDr;auxDY:=auxDY+intT*PsiDi; end
         else begin { ********** en-dessous du dioptre  ********************  }
              end
         end;
      deuxd.Spect[indModU].data[i,j,0].x:=X/Xunit;deuxd.Spect[indModU].data[i,j,0].y:=Y/Yunit;
      deuxd.Spect[indModU].data[i,j,0].z:=sqr(auxUX)+sqr(auxUY);
      deuxd.Spect[indModD].data[i,j,0].x:=X/Xunit;deuxd.Spect[indModD].data[i,j,0].y:=Y/Yunit;
      deuxd.Spect[indModD].data[i,j,0].z:=sqr(auxDX)+sqr(auxDY);
      if PhaseYN then begin
      deuxd.Spect[indArgU].data[i,j,0].x:=X/Xunit;deuxd.Spect[indArgU].data[i,j,0].y:=Y/Yunit;
      deuxd.Spect[indArgU].data[i,j,0].z:=Arg(Cplx(auxUX,auxUY))/rad; end;
    end;
    Application.processMessages;if Stop then break;
    ProgressBar.position:=i;ProgressBar.update;
  end;
  status.Caption:='Idle';status.Update;
end;


procedure TWaveFunctionPlotting.PolarisedNeutronsWithSFbuttonClick(Sender: TObject);
begin
   initParameters;generateDistributions;
   PolarisedNeutronsWithSF('Phi');
   sendmessage(principalform.handle,wm_redessine,0,0);
end;

procedure TWaveFunctionPlotting.CutOnDomainsUpDown1D(nameD,nameP:string);
var i,m,n:integer;
    NbT:integer;
    indMod,indArg,indD,indCtrl:integer;
    X,Y,Xs,PsiX,PsiY,pol,auxX,auxY,auxI:double;
    p,q,q1,k1,R,T,Rx,Ry,Xlength,intT:double;
    phaseYN:boolean;
begin
   indD:=deuxd.indexof(nameD);
   indMod:=und.addspectre(nameP+'_mod',CutSize,1,1);
   PhaseYN:=PhaseBox.checked;
   if PhaseYN then begin indArg:=und.addspectre(nameP+'_arg',CutSize,1,1);
                         indCtrl:=und.addspectre('Ctrl',CutSize,1,1); end;
   Xlength:=Xmax-Xmin;
   NbT:=trunc(NbPtsT/2);
   for i:=0 to CutSize-1 do
   begin
       X:=Xc;Y:=i*Ystep;
       auxX:=0;auxY:=0;
       for m:=-NbT to NbT do
         begin
         if (Y>0)
         then begin
            intT:=DistrTint[m+NbT];
            p:=k0*DistrTc[m+NbT]; { composante k //    }
            q:=k0*DistrTs[m+NbT]; { composante k perp  }
            { d'ou vient la réflexion ? }
            Xs:=X-p/q*Y; { abscisse de l'origine du rayon (pour Y0=0) }
{            While (Xs<0)    do Xs:=Xs+Xlength; { cela ne donne pas de très jolis résultats }
{            While (Xs>Xmax) do Xs:=Xs-Xlength; }
            if (Xs<Xmin) or (Xs>Xmax) { inutile }
            then begin { onde plane seulement }
                 auxI:=p*(X-X0)-q*(Y-Y0);
                 auxX:=auxX+intT*cos(auxI);auxY:=auxY+intT*sin(auxI);
                 end
             else begin
                 { quelle est la polarisation du domaine reflecteur ? }
                 pol:=deuxd.spect[indD].data[trunc((Xs-Xmin)/Xstep),0,0].z;
                 if pol=1 then k1:=indiceU*k0 else k1:=indiceD*k0;
                 calculPsiR(X,Y,p,q,k1,0,PsiX,PsiY);
                 auxX:=auxX+intT*PsiX;auxY:=auxY+intT*PsiY;
                 end;
            end;
         end;
       und.Spect[indMod].data[i,0,0].x:=Y/Yunit;
       und.Spect[indMod].data[i,0,0].z:=sqr(auxX)+sqr(auxY); {Cabs(aux); }
       if PhaseYN then begin
       und.Spect[indArg].data[i,0,0].x:=Y/Yunit;
       und.Spect[indArg].data[i,0,0].z:=Arg(Cplx(auxX,auxY))/rad; { repassage en deg }
       und.Spect[indCtrl].data[i,0,0].x:=Y/Yunit;
       und.Spect[indCtrl].data[i,0,0].z:=pol; end;
   end;
end;

procedure TWaveFunctionPlotting.CutafterreflectiononmagneticdomainsUPDO1Click(Sender: TObject);
var i,indMod:integer;
begin
   initParameters;generateDistributions;
   GenerateUPDOWNpattern('DomainsUD');
   CutOnDomainsUpDown1D('DomainsUD','Phi');
   for i:=2 to NrDomainsAveragesBox.value do
   begin
     GenerateUPDOWNpattern('DomainsUD');
     CutOnDomainsUpDown1D('DomainsUD','Temp');
     und.Addition('Phi_mod','Temp_mod','Phi_mod')
   end;
   sendmessage(principalform.handle,wm_redessine,0,0);
end;

{ test de la reflexion sur des domaines }
procedure TWaveFunctionPlotting.OnDomainsUPDOWN(nameD,nameP:string);
var i,j,m,n:integer;
    NbT:integer;
    indMod,indArg,indD,indCtrl:integer;
    X,Y,Xs,PsiX,PsiY,pol,auxX,auxY,auxI:double;
    p,q,q1,k1,R,T,Rx,Ry,Xlength,intT:double;
    phaseYN:boolean;
begin
   stop:=false;
   status.Caption:='Calculating...';status.Update;
   ProgressBar.max:=Xsize;ProgressBar.position:=0;

   indD:=deuxd.indexof(nameD);
   indMod:=deuxd.addspectre(nameP+'_mod',Xsize,Ysize,1);
   PhaseYN:=PhaseBox.checked;
   if PhaseYN then begin indArg:=deuxd.addspectre(nameP+'_arg',Xsize,Ysize,1);
                         indCtrl:=deuxd.addspectre('Ctrl',Xsize,Ysize,1); end;

   Xlength:=Xmax-Xmin;
   NbT:=trunc(NbPtsT/2);
   for i:=0 to Xsize-1 do
   begin
     ProgressBar.position:=i;ProgressBar.update;
     for j:=0 to Ysize-1 do
     begin
       X:=xmin+i*Xstep;Y:=ymin+j*Ystep;
       auxX:=0;auxY:=0;
       for m:=-NbT to NbT do
         begin
         if (Y>0)
         then begin
            intT:=DistrTint[m+NbT];
            p:=k0*DistrTc[m+NbT]; { composante k //    }
            q:=k0*DistrTs[m+NbT]; { composante k perp  }
            { d'ou vient la réflexion ? }
            Xs:=X-p/q*Y; { abscisse de l'origine du rayon (pour Y0=0) }
            While (Xs<0)    do Xs:=Xs+Xlength; { cela ne donne pas de très jolis résultats }
            While (Xs>Xmax) do Xs:=Xs-Xlength;
            if (Xs<Xmin) or (Xs>Xmax) { inutile }
            then begin { onde plane seulement }
                 auxI:=p*(X-X0)-q*(Y-Y0);
                 auxX:=auxX+intT*cos(auxI);auxY:=auxY+intT*sin(auxI);
                 end
            else begin
                 { quelle est la polarisation du domaine reflecteur ? }
                 pol:=deuxd.spect[indD].data[trunc((Xs-Xmin)/Xstep),0,0].z;
                 if pol=1 then k1:=indiceU*k0 else k1:=indiceD*k0;
                 calculPsiR(X,Y,p,q,k1,0,PsiX,PsiY);
                 auxX:=auxX+intT*PsiX;auxY:=auxY+intT*PsiY;
                 end;
            end;
         end;
       deuxd.Spect[indMod].data[i,j,0].x:=X/Xunit;deuxd.Spect[indMod].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indMod].data[i,j,0].z:=sqr(auxX)+sqr(auxY); {Cabs(aux); }
       if PhaseYN then begin
       deuxd.Spect[indArg].data[i,j,0].x:=X/Xunit;deuxd.Spect[indArg].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indArg].data[i,j,0].z:=Arg(Cplx(auxX,auxY))/rad; { repassage en deg }
       deuxd.Spect[indCtrl].data[i,j,0].x:=X/Xunit;deuxd.Spect[indCtrl].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indCtrl].data[i,j,0].z:=pol; end;
     end;
     Application.processMessages;if Stop then break;
   end;
   status.Caption:='Idle';status.Update;
end;

{ test de la reflexion sur des domaines }
procedure TWaveFunctionPlotting.OnDomainesUPDOWNClick(Sender: TObject);
begin
   initParameters;generateDistributions;
   GenerateUPDOWNpattern('DomainsUD');
   OnDomainsUpDown('DomainsUD','Phi');
   sendmessage(principalform.handle,wm_redessine,0,0);
end;

{ A FAIRE }
procedure TWaveFunctionPlotting.OnRandomDomainsButtonClick(Sender: TObject);
begin
 { clkqfljdls }
end;

{ simulation d'une série de domaines UP et DOWN                                 }
{ on fait l'hypothèse que l'aimantation est parallèle au champ magnétique appliqué }
{ on ne génère que des domaines UP (+1) et (DOWN)                               }
procedure TWaveFunctionPlotting.GenerateUPDOWNpattern(name:string);
var i,j,iCrt,iMax,ind:integer;
    S,pol:integer; { polarisation des domaines }
    X,Xd,Y,sizeD,dSizeD:double;
begin
  { création de la figure des domains }
  ind:=deuxd.addspectre(name,Xsize,Ysize,1);
  SizeD:=sizeDbox.value;dSizeD:=dSizeDBox.value;
  iCrt:=0;pol:=-1;X:=xmin;
  while iCrt<Xsize do
  begin
    S:=round(randG(SizeD,dSizeD));
    Xd:=X+S;
    if pol=1 then pol:=-1 else pol:=1;
    while (X<Xd) and (iCrt<Xsize) do
    begin
      for j:=0 to Ysize-1 do
        begin
        Y:=ymin+j*Ystep;
        deuxd.spect[ind].data[iCrt,j,0].x:=X/Xunit;
        deuxd.spect[ind].data[iCrt,j,0].y:=Y/Yunit;
        if (Y>0)  then deuxd.spect[ind].data[iCrt,j,0].z:=0;
        if (Y<=0) then deuxd.spect[ind].data[iCrt,j,0].z:=pol;
        end;
      iCrt:=iCrt+1;
      X:=xmin+iCrt*Xstep;
    end;
  end;
end;

{ simulation d'une série de domaines UP et DOWN                                 }
{ on fait l'hypothèse que l'aimantation est parallèle au champ magnétique appliqué }
{ on ne génère que des domaines UP (+1) et (DOWN)                               }
procedure TWaveFunctionPlotting.GenerateUPDOWNpatternButtonClick(Sender: TObject);
begin
  initParameters;generateDistributions;
  GenerateUPDOWNpattern('DomainsUD');
  sendmessage(principalform.handle,wm_redessine,0,0);
end;

procedure TWaveFunctionPlotting.GeneraterandomPattern(name:string);
var i,j,iCrt,iMax,ind:integer;
    S:integer; { polarisation des domaines }
    Phi,phi0,dPhi,X,Xd,Y,sizeD,dSizeD:double;
begin
  { création de la figure des domains }
  ind:=deuxd.addspectre('Domains',Xsize,Ysize,1);
  SizeD:=sizeDbox.value;dSizeD:=dSizeDBox.value;
  Phi0:=phibox.value;dPhi:=dPhiBox.value;
  iCrt:=0;X:=xmin;
  while iCrt<Xsize do
  begin
    S:=round(randG(SizeD,dSizeD));
    phi:=round(randG(Phi0,dPhi)+36000) mod 360; { angle de l'aimantation }
    Xd:=X+S;
    while (X<Xd) and (iCrt<Xsize) do
    begin
      for j:=0 to Ysize-1 do
        begin
        Y:=ymin+j*Ystep;
        deuxd.spect[ind].data[iCrt,j,0].x:=X/Xunit;
        deuxd.spect[ind].data[iCrt,j,0].y:=Y/Yunit;
        if (Y>0)  then deuxd.spect[ind].data[iCrt,j,0].z:=0;
        if (Y<=0) then deuxd.spect[ind].data[iCrt,j,0].z:=phi;
        end;
      iCrt:=iCrt+1;
      X:=xmin+iCrt*Xstep;
    end;
  end;
end;

procedure TWaveFunctionPlotting.GeneraterandomPatternButtonClick(Sender: TObject);
begin
  initParameters;generateDistributions;
  GenerateRandomPattern('Domains');
  sendmessage(principalform.handle,wm_redessine,0,0);
end;

{ ***************************************************************************** }
{ ***************************************************************************** }
{   FINITE SIZE                                                                 }
{ ***************************************************************************** }
{ ***************************************************************************** }

{ calcul de la réflexion d'une onde plane sur un objet de taille finie          }
procedure TWaveFunctionPlotting.Planewave1Click(Sender: TObject);
var i,j,m,n:integer;
    NbT:integer;
    indMod,indArg:integer;
    X,Y,auxI,auxR,auxT,auxX,auxY,aux,aux1:double;
    p,q,q1,k1,R,T,Rx,Ry:double;
    ThetaIn,TTin,L,length:double; { longueur de l'objet }
begin
   stop:=false;
   status.Caption:='Calculating...';status.Update;
   initParameters;generateDistributions;
   length:=lengthBox.value; { longueur de l'objet }
   L:=length/2;
   ProgressBar.max:=Xsize;ProgressBar.position:=0;

   indMod:=deuxd.addspectre('Phi_mod',Xsize,Ysize,1);
   indArg:=deuxd.addspectre('Phi_arg',Xsize,Ysize,1);

   { on ne tient pas compte d'une quelconque distribution }
   ThetaIn:=Theta0box.value*pi/180; { angle d'incidence }
   TTin:=tan(thetaIn);
   for i:=0 to Xsize-1 do
   begin
     ProgressBar.position:=i;ProgressBar.update;
     for j:=0 to Ysize-1 do
     begin
       X:=xmin+i*Xstep;Y:=ymin+j*Ystep;
       auxX:=0;auxY:=0;
       { calcul de la fonction d'onde }
       p:=k0*cos(thetaIn); { composante k //    }
       q:=-k0*sin(thetaIn); { composante k perp  }
       k1:=indice*k0;        { vecteur d'onde dans le dioptre     }
       { est-ce qu'on est dans la zone ou une onde est réfléchie ? }
       if (Y>0) and (Y<((X+L)*TTin)) and (Y>((X-L)*TTin))
       then begin
            if (k1>=p)
            then begin { réflexion partielle } { ok }
                 q1:=sqrt(sqr(k1)-sqr(p)); { vecteur d'onde réel }
                 R:=(abs(q)-q1)/(abs(q)+q1);  { coefficient de réflexion réel ; pas de facteur de phase...curieux... }
                 auxI:=p*(X-X0)+q*(Y-Y0);             { facteur de phase de l'onde incidente }
                 auxR:=p*(X-X0)-q*(Y-Y0);             { facteur de phase de l'onde réfléchie }
                 auxX:=cos(auxI)+R*cos(auxR);
                 auxY:=sin(auxI)+R*sin(auxR); end
            else begin { onde totalement réfléchie } {ok }
                 q1:=sqrt(-sqr(k1)+sqr(p)); { composante du vecteur d'onde imaginaire pure }
                 aux:=(sqr(q)+sqr(q1));
                 Rx:=(sqr(q)-sqr(q1))/aux;Ry:=-2*q*q1/aux; { CHECK SIGNS }
                 auxI:=p*(X-X0)+q*(Y-Y0);             { facteur de phase de l'onde incidente }
                 auxR:=p*(X-X0)-q*(Y-Y0);
                 auxX:=cos(auxI)+Rx*cos(auxR)-Ry*sin(auxR);
                 auxY:=sin(auxI)+Rx*sin(auxR)+Ry*cos(auxR); end;
            end
       else begin { on est dans la zone ou il y a simplement une onde plane  }
            auxR:=k0*(cos(thetaIn)*(X-X0)-sin(thetaIn)*(Y-Y0));
            auxX:=cos(auxR);
            auxY:=sin(auxR);
            end;

       deuxd.Spect[indMod].data[i,j,0].x:=X/Xunit;
       deuxd.Spect[indMod].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indMod].data[i,j,0].z:=sqr(auxX)+sqr(auxY); {Cabs(aux); }
       deuxd.Spect[indArg].data[i,j,0].x:=X/Xunit;
       deuxd.Spect[indArg].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indArg].data[i,j,0].z:=Arg(Cplx(auxX,auxY))/rad; { repassage en deg }
     end;
     Application.processMessages;if Stop then break;
   end;
   status.Caption:='Idle';status.Update;
end;

{ calcul de la reflection sur n objet de taille finie                           }
procedure TWaveFunctionPlotting.FiniteSizeTransversecoherenceOnlyClick(Sender: TObject);
var i,j,m,n:integer;
    NbT:integer;
    indMod,indArg:integer;
    X,Y,auxI,auxR,auxT,auxX,auxY,aux,aux1:double;
    p,q,q1,k1,R,T,Rx,Ry:double;
    length,L,TTin,TTrefracted:double; { longueur de l'objet }
begin
   stop:=false;
   status.Caption:='Calculating...';status.Update;
   initParameters;generateDistributions;
   length:=lengthBox.value; { longueur de l'objet }
   L:=length/2;
   ProgressBar.max:=Xsize;ProgressBar.position:=0;

   indMod:=deuxd.addspectre('Phi_mod',Xsize,Ysize,1);
   indArg:=deuxd.addspectre('Phi_arg',Xsize,Ysize,1);

   NbT:=trunc(NbPtsT/2);
   for i:=0 to Xsize-1 do
   begin
     ProgressBar.position:=i;ProgressBar.update;
     for j:=0 to Ysize-1 do
     begin
       X:=xmin+i*Xstep;Y:=ymin+j*Ystep;
       auxX:=0;auxY:=0;
       for m:=-NbT to NbT do
         begin
            { calcul de la fonction d'onde }
            p:=k0*DistrTc[m+NbT]; { composante k //    }
            q:=k0*DistrTs[m+NbT]; { composante k perp  } { QUESTION SUR LE SIGNE }
            TTin:=q/p; { tan(theta in) }
            k1:=indice*k0;        { vecteur d'onde dans le dioptre     }
            if (k1>=p) { reflexion partielle }
            then begin q1:=sqrt(sqr(k1)-sqr(p));TTrefracted:=q1/p end
            else TTrefracted:=0; { relexion totale }
            if (Y>0) { au-dessus de la couche }
            then if (Y<((X+L)*TTin)) and (Y>((X-L)*TTin))
                 then begin { ********** au-dessus du dioptre  ********************  }
                      if (k1>=p)
                      then begin { réflexion partielle } { ok }
                      q1:=sqrt(sqr(k1)-sqr(p)); { vecteur d'onde réel }
                      R:=(q-q1)/(q+q1);  { coefficient de réflexion réel ; pas de facteur de phase...curieux... }
                      auxI:=p*(X-X0)-q*(Y-Y0);             { facteur de phase de l'onde incidente }
                      auxR:=p*(X-X0)+q*(Y-Y0);             { facteur de phase de l'onde réfléchie }
                      auxX:=auxX+DistrTint[m+NbT]*(cos(auxI)+R*cos(auxR));
                      auxY:=auxY+DistrTint[m+NbT]*(sin(auxI)+R*sin(auxR)); end
                      else begin { onde totalement réfléchie } {ok }
                      q1:=sqrt(-sqr(k1)+sqr(p)); { composante du vecteur d'onde imaginaire pure }
                      aux:=(sqr(q)+sqr(q1));
                      Rx:=(sqr(q)-sqr(q1))/aux;Ry:=-2*q*q1/aux;
                      auxI:=p*(X-X0)-q*(Y-Y0);             { facteur de phase de l'onde incidente }
                      auxR:=p*(X-X0)+q*(Y-Y0);
                      auxX:=auxX+DistrTint[m+NbT]*(cos(auxI)+Rx*cos(auxR)-Ry*sin(auxR));
                      auxY:=auxY+DistrTint[m+NbT]*(sin(auxI)+Rx*sin(auxR)+Ry*cos(auxR)); end
                      end
                 else begin { onde plane seulement }
                      auxI:=p*(X-X0)-q*(Y-Y0);             { facteur de phase de l'onde incidente }
                      auxX:=auxX+DistrTint[m+NbT]*cos(auxI);
                      auxY:=auxY+DistrTint[m+NbT]*sin(auxI);
                      end
            { sous la surface ; le vecteur d'onde peut etre evanescent !!! }
            else begin
                 if (Y>(-(X+L)*TTrefracted)) and (Y<(-(X-L)*TTrefracted))
                 then begin { ********** en-dessous du dioptre  ********************  }
                      aux:=DistrTint[m+NbT];
                      { calcul du vecteur d'onde dans le milieu   }
                      if (k1>=p)
                      then begin { réflexion partielle } { ok }
                           q1:=sqrt(+sqr(k1)-sqr(p));
                           T:=DistrTint[m+NbT]*2*q/(q+q1); { coefficient de transmission }
                           auxT:=p*(X-X0)-q1*(Y-Y0);
                           auxX:=auxX+T*cos(auxT);
                           auxY:=auxY+T*sin(auxT); end
                      else begin { onde evanescente } { ok }
                           q1:=sqrt(-sqr(k1)+sqr(p)); { composantes imaginaire pure ! }
                           aux:=DistrTint[m+NbT]*2*q/(sqr(q)+sqr(q1))*exp(q1*y);
                           auxT:=p*(X-X0);
                           auxX:=auxX+aux*(q*cos(auxT)+q1*sin(auxT));
                           auxY:=auxY+aux*(q*sin(auxT)-q1*cos(auxT)); end;
                      end;
                 if (Y<(-(X+L)*TTin)) or (Y>(-(X-L)*TTin))
                 then begin { onde plane seulement }
                      auxI:=p*(X-X0)-q*(Y-Y0);             { facteur de phase de l'onde incidente }
                      auxX:=auxX+DistrTint[m+NbT]*cos(auxI);
                      auxY:=auxY+DistrTint[m+NbT]*sin(auxI);
                      end
                 end;
         end;
       deuxd.Spect[indMod].data[i,j,0].x:=X/Xunit;
       deuxd.Spect[indMod].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indMod].data[i,j,0].z:=sqr(auxX)+sqr(auxY); {Cabs(aux); }
       deuxd.Spect[indArg].data[i,j,0].x:=X/Xunit;
       deuxd.Spect[indArg].data[i,j,0].y:=Y/Yunit;
       deuxd.Spect[indArg].data[i,j,0].z:=Arg(Cplx(auxX,auxY))/rad; { repassage en deg }
     end;
     Application.processMessages;if Stop then break;
   end;
   status.Caption:='Idle';status.Update;
end;


{ ***************************************************************************** }
{ ***************************************************************************** }
{ MOVIES }
{ ***************************************************************************** }
{ ***************************************************************************** }

procedure TWaveFunctionPlotting.Save1Click(Sender: TObject);
begin
  // Set file extension and filter to GIF (or whatever is currently displayed)
  SaveDialog.DefaultExt := GraphicExtension(TGraphicClass(ImageAnimate.Picture.Graphic.ClassType));
  SaveDialog.Filter :=  GraphicFilter(TGraphicClass(ImageAnimate.Picture.Graphic.ClassType));
  // Prompt for filename
  if (SaveDialog.Execute) then
    ImageAnimate.Picture.SaveToFile(SaveDialog.Filename);
end;

procedure TWaveFunctionPlotting.Animate1Click(Sender: TObject);
var
  GIF			: TGIFImage;
  i, Tsize, indMod     	: integer;
  OptimizeOptions	: TGIFOptimizeOptions;
  tmin,tmax,tstep:double;
  screenCopy:Tbitmap;
  crect : Trect;
  nom:string;

  function TransparentIndex(GIF: TGIFSubImage): byte;
  begin
    // Use the lower left pixel as the transparent color
    Result := GIF.Pixels[0, GIF.Height-1];
  end;

  function AddBitmap(GIF: TGIFImage; Source: TGraphic; Transparent: boolean): integer;
  var
    Ext			: TGIFGraphicControlExtension;
    LoopExt		: TGIFAppExtNSLoop;
  begin
    // Preview the image being added (for user feedback)
    ImageAnimate.Picture.Assign(Source);
    ImageAnimate.Update;

    // Add the source image to the animation
    Result := GIF.Add(Source);
    // Netscape Loop extension must be the first extension in the first frame!
    if (Result = 0) then
    begin
      LoopExt := TGIFAppExtNSLoop.Create(GIF.Images[Result]);
      LoopExt.Loops := 0; // Number of loops (0 = forever)
      GIF.Images[Result].Extensions.Add(LoopExt);
    end;
    // Add Graphic Control Extension
    Ext := TGIFGraphicControlExtension.Create(GIF.Images[Result]);
    Ext.Delay := 30; // Animation delay (30 = 300 mS)
    if (Transparent) then
    begin
      Ext.Transparent := True;
      Ext.TransparentColorIndex := TransparentIndex(GIF.Images[Result]);
    end;
    GIF.Images[Result].Extensions.Add(Ext);
  end;

begin
  { variables temporelles }
  tmin:=TminBox.value;tstep:=TstepBox.value;tmax:=TmaxBox.value;Tsize:=round((Tmax-Tmin)/Tstep)+1;
  try
    GIF := TGIFImage.Create;
    try
      { balayage des courbes calculées }
      for i := 1 to Tsize do
      begin
        nom:='Phi'+IntToStr(i)+'_mod';
        Graph2D3DForm.ChoixSpectre.Text:=nom;
        Graph2D3DForm.LoadSpectrum(nom);
        Graph2D3DFOrm.Refresh;
        Graph2D3Dform.copyToClipboard;
        Image1.Picture.Assign(ClipBoard);
{        GIF.add(Image1.Picture.Graphic); }
        AddBitmap(GIF, Image1.picture.graphic, false);
      end;
      // Optimize Color map...
      if (CheckBoxPalette.Checked) then GIF.OptimizeColorMap;
      // Optimize GIF frames...
      OptimizeOptions := [];
      if (CheckBoxMerge.Checked) then include(OptimizeOptions, ooMerge);
      if (CheckBoxCrop.Checked)  then include(OptimizeOptions, ooCrop);
      if (OptimizeOptions <> []) then GIF.Optimize(OptimizeOptions, rmNone, dmNearest, 0);
      // Display the animation
      ImageAnimate.Picture.Assign(GIF);
    finally
      GIF.Free;
    end;
  finally
    Screen.Cursor := crDefault;
  end;
end;

procedure TWaveFunctionPlotting.Anime1Click(Sender: TObject);
var i,ind,Tisze:integer;
    tmin,tstep,tmax:double;
    nom:string;
    temp:Tbitmap;
begin
    temp:=TBitmap.create;
    tmin:=TminBox.value;tstep:=TstepBox.value;tmax:=TmaxBox.value;Tsize:=round((Tmax-Tmin)/Tstep)+1;
    for i := 1 to Tsize do
    begin
        nom:='Phi'+IntToStr(i)+'_mod';
        Graph2D3DForm.ChoixSpectre.Text:=nom;
        Graph2D3DForm.LoadSpectrum(nom);
        Graph2D3DFOrm.Refresh;
        Graph2D3Dform.copyToClipboard;
        Image1.Picture.Assign(ClipBoard);
        ind:=ImageList.add(Image1.Picture.Bitmap,nil);
        ImageList.GetBitmap(ind,temp); { marche pas ... }
        ImageAnimate.Picture.Assign(temp);
        showMessage(IntToStr(ind));
    end;
end;

procedure TWaveFunctionPlotting.Graph2D3DFormStartMenu1Click(Sender: TObject);
begin
  JvBmpAnimator.Active:=true;
end;



end.


