% File: mypde_3_main1_with_PhononScattering_ElasticScattering_radioactive_and_augerdecay.m

%clear all
clc
tic
% Global area
global nr  ne  dr de r e  N el rl N0 c  M c1 c2 c6  ncall
%
% Model parameters

el=100;      
rl=200;       

Mx=2.61*9.109*10^(-31);  %Exciton mass in kg
%a=0.5 micro ev/micro_m^2
M=10;       
v=4500; %Speed of sound in m/s  
h1=1.055*10^(-34);  %Reduced Planck constant in Joul.second 
c3=(4*(Mx^2)*(v^2))/((h1^2)*(M));
e0=601.04;           %micro_ev
a=0.801*10^(-25);  % Steepness constant in unit Joul/micro_m^2
kB=1.381*10^(-23);  % Boltzmann constant in unit J/k 
alp=(0.5*M)/(e0);

c1=(9)/(sqrt(M));          %c1=((2*v)/sqrt(M)) [unit=micro_m/ns]
c2=0.01497*(sqrt(M)); %c2=(4*(a/e0)*v*sqrt(M)) [unit=1/(micro_m*ns)]

c4=(1/((2*pi)^3))*((((4*pi)^2)*sqrt(2)*(Mx^(3/2)))/(h1^3))*(((e0*1.602*10^(-25))/M)^(3/2))*(10^(-18));

h=6.972;            %h=e0/K_B=(2*Mx*v^2)/k_B  [unit=kelvin]
tem=3;        % Temperature in Kelvin
c=(h/(tem*M));

n0=5.85*10^(-4);        


c6=95.275;        %[unit=1/ns][c6=(4*((al)^2)*Mx*(e0^2))/(pi*(h1^3))] [a1=scattering length=2 nm, Mx=exciton mass=2.61*9.109*10^(-31)kg, e0=0.6 meV,pi=3.14, h1=Reduced Planck constant=1.055*10^(-34)]
fprintf('check1');
% Grid in k direction
  ne=100;           
  de=el/(ne-1);
  %de=((e0*10^(-3))/M);
  for j=1:ne
    e(j)=(j-1)*de;
  end
%
% Grid in r direction
  nr=100;              
  dr=rl/(nr-1);
  %dr=2;
  for i=1:nr
    r(i)=(i-1)*dr;
  end
 fprintf('check2'); 
%
% Independent variable time for ODE integration
  tf=200;
  tout=(0:1:tf)';
  nout=201;
  ncall=0;
  
  
% Preallocation for speed
N0=zeros(nr,ne);
y0=zeros(nr,ne);
N_rad_1=zeros(nr,ne);
N_rad_5=zeros(nr,ne);
N_rad_10=zeros(nr,ne);
N_rad_50=zeros(nr,ne);
N_rad_100=zeros(nr,ne);
N_rad_150=zeros(nr,ne);
N_rad_200=zeros(nr,ne);
fprintf('check3'); 

% Initial condition
  for i=1:nr
  for j=1:ne
      if (i>=35) 
          N0(i,j)=0;
      else
          N0(i,j)=n0*(exp((-(j-(45-(alp.*((r(i)).^2)))).^2)/(2*(2)^2)));
          
      end 
          y0((i-1)*ne+j)=N0(i,j);
      
  end
  end

%
% ODE integration 
  
  
  [t,y]=ode15s(@programme_BTE,tout,y0);    
 
  fprintf('check4'); 
  
  
%
%   2D to 3D matrices
  for it=1:nout
  for i=1:nr
  for j=1:ne 
    N(it,i,j)=y(it,(i-1)*ne+j);
    N(1,i,j)=y0((i-1)*ne+j);
    N(1,i,j)=N0(i,j);
    
  end
  end
  end
 
fprintf('check5');   
  

 
% Radial profiles (at t = tf)

for i=1:nr
  for j=1:ne
    N_rad_1(i,j)=N(1,i,j);
  end
end
 
  for i=1:nr
  for j=1:ne
    N_rad_5(i,j)=N(5,i,j);
  end
  end
  
   for i=1:nr
  for j=1:ne
    N_rad_10(i,j)=N(10,i,j);
  end
  end
  
  for i=1:nr
  for j=1:ne
    N_rad_50(i,j)=N(50,i,j);
  end
  end
  
  for i=1:nr
  for j=1:ne
    N_rad_100(i,j)=N(100,i,j);
  end
  end
  
  for i=1:nr
  for j=1:ne
    N_rad_150(i,j)=N(150,i,j);
  end
  end
  
  for i=1:nr
  for j=1:ne
    N_rad_200(i,j)=N(200,i,j);
  end
  end
  
  % Summation of exciton number for a particular time


  N_sum_1=0;
  for i=1:nr
      for j=1:ne
          N_sum_1=N_sum_1+(c4*((r(i))^2).*sqrt(j).*N_rad_1(i,j).*dr);
          
      end
  end 
 fprintf('\n N_sum_1 = %5d\n',N_sum_1);
  
 
 
  N_sum_5=0;
  for i=1:nr
      for j=1:ne
          N_sum_5=N_sum_5+(c4*((r(i))^2).*sqrt(j).*N_rad_5(i,j).*dr);
          
      end
  end 
  fprintf('\n N_sum_5 = %5d\n',N_sum_5);
  
   N_sum_10=0;
  for i=1:nr
      for j=1:ne
          N_sum_10=N_sum_10+(c4*((r(i))^2).*sqrt(j).*N_rad_10(i,j).*dr);
          
      end
  end 
  fprintf('\n N_sum_10 = %5d\n',N_sum_10);
  
   N_sum_50=0;
  for i=1:nr
      for j=1:ne
          N_sum_50=N_sum_50+(c4*((r(i))^2).*sqrt(j).*N_rad_50(i,j).*dr);
          
      end
  end 
 fprintf('\n N_sum_50 = %5d\n',N_sum_50);
  
   N_sum_100=0;
  for i=1:nr
      for j=1:ne
          N_sum_100=N_sum_100+(c4*((r(i))^2).*sqrt(j).*N_rad_100(i,j).*dr);
          
      end
  end 
  fprintf('\n N_sum_100 = %5d\n',N_sum_100);
  
   N_sum_150=0;
  for i=1:nr
      for j=1:ne
          N_sum_150=N_sum_150+(c4*((r(i))^2).*sqrt(j).*N_rad_150(i,j).*dr);
          
      end
  end 
  fprintf('\n N_sum_150 = %5d\n',N_sum_150);
  
  
    N_sum_200=0;
  for i=1:nr
      for j=1:ne
         N_sum_200=N_sum_200+(c4*((r(i))^2).*sqrt(j).*N_rad_200(i,j).*dr);
          
      end
  end 
  fprintf('\n N_sum_200 = %5d\n',N_sum_200);
 
  


% 3D plotting
  
  
  figure()
  surf(e,r,N_rad_1)
  figure()
  surf(e,r,N_rad_5)
  
  figure()
  surf(e,r,N_rad_10)
  figure()
  surf(e,r,N_rad_50)
  figure()
  surf(e,r,N_rad_100)
  figure()
  surf(e,r,N_rad_150)
  
  figure()
  surf(e,r,N_rad_200)
 
  toc
  save 'C:\MatlabProgram\trap_ps_es_condensate_tbd_3K_M10_n0_5point85times10^(-4)_200ns_100points.mat'