%change directory to this file folder: cd('')
clear;
clc;
path = 'D:\OneDrive - gznu.edu.cn\Project\cee\anode_strip\Prob3\';
tic;
pwd
import com.comsol.model.*
import com.comsol.model.util.*
modelName = 'BM_ro_strip';
model = mphload(modelName);

% mobility  300V/cm,   4.5cm/us
mu = 4.5/300;   % unit  cm^2/v/us
% x0, y0, z0:    The start point of one electron,   unit mm
x0 = 0;
y0 = 0;
z0 =3.5;
step = 0.01; % step length 0.5mm
% the position of the strip plan
pcb_t = 1.234;
strip_t = 0.01;
z_ref = pcb_t+2*strip_t;     % the height of the strip plan    


% begin the drift 
i = 1; %  the ith step
x(i) = x0; y(i) = y0; z(i) = z0;


while(1)
%     fprintf('i = %d\n',i);
  %  extract Electric Field from COMSOL 
%     fprintf('x: %g\t y: %g\t  z: %g\n', x(i), y(i), z(i));
    if z(i)<z_ref
        % collected by electrode
        break;
    end
    % velocity of electron  v = mu * E,  E:V/cm, v:cm/microS
    [Ex,Ey,Ez,unit]= mphinterp(model,{'-es.Ex','-es.Ey','-es.Ez'},'coord',[x(i),y(i),z(i)]');% V/m
    % flee from the cell 
	if isnan(Ex)|| isnan(Ey)|| isnan(Ez)
        break;
    end
    Vx = mu*Ex *1.0e-1; Vy = mu*Ey *1.0e-1;   Vz = mu*Ez *1.0e-1;          % V: mm/microS
%     fprintf('Ex: %g\t Ey: %g\t Ez: %g \n', Ex, Ey, Ez);
%     fprintf('Vx: %g\t Vy: %g\t Vz: %g\n', Vx, Vy, Vz);

    % weighting field
    [Ewx,Ewy,Ewz,unit]= mphinterp(model,{'es3.Ex','es3.Ey','es3.Ez'},'coord',[x(i),y(i),z(i)]');% V/m
%      fprintf('Ewx: %g\t Ewy: %g\t Ewz: %g \n', Ewx, Ewy, Ewz);
    % induced current of electron q = -1
    current(i) = - (Ewx*Vx + Ewy*Vy+ Ewz*Vz)*1.0e3; % A
    % time step 
    dt(i) = step/sqrt(Vx^2+ Vy^2+ Vz^2); %us
    % fprintf('dt: %g\n', dt(i));
    charge(i) = current(i)*dt(i)*1.0e-6; % C      
    x(i+1) = x(i) + (dt(i)*Vx); y(i+1) = y(i) + dt(i)*Vy; z(i+1) = z(i) + dt(i)*Vz;% mm
    i = i + 1;
end

% charge vs z
format = '.fig';
charge_y_name = 'induced_current VS z';
charge_y_name = [path charge_y_name format];
charge_y_plot = figure(1);
set(charge_y_plot, 'Position', [10 10 1000 400])
subplot(1,2,1);
plot(z(1:length(z)-1), current,'*');
xlabel('z [mm]');
ylabel('induced current of one electron[A]');
set(gca,'xdir','reverse');
grid on;
subplot(1,2,2);
plot(cumsum(dt*1.0e3), current,'*');
charge_sum = sum(current.*dt*1e-6);

xlabel(['the drift time [ns]', '     charge = ', num2str(charge_sum),' [C]'])
ylabel('induced current of one electron [A]');

grid on;
saveas(charge_y_plot, charge_y_name);


model.save(modelName);

toc











































