% Criating file with n, V, T, Q, Pot and E profiles, and plotting them


clear all
close all

pp= 'b_B.';               % Profiles <name> without .inp*

%%%%%%%%%%%%%% Input %%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ntot= 48;               % Total number of Procces.
Lx= 0.001;              % System size [m]
mm= 0;                  % Flag for the magnetic field (if = 0 - no magnetic field)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pn= 'inp.dat.';
pn= strcat(pp,pn);
vout=  strcat('pr_',pp,'txt');

%%%%%%%%%%% Criate variable %%%%%%%%%%%%%%%%%%%%%

N=Ntot-1;
N0= 0;

for i=N0:N
    pp=strcat(pn,num2str(i))
    par = load(pp); 
    if (i==N0)
        sL= size(par,1);
        sr= size(par,2);
        fsz= (sL-1).*(N-N0 + 1)+1;
        de= zeros(fsz,sr);
        for j= 1:sL;
            de(j,:)= par(j,:);
        end
    else
        zz= (sL-1).*(i-N0)+1;
        par(1,sr-1)= 0.0;
        par(1,sr)= 0.0;
        de(zz,:)= de(zz,:) + par(1,:);
           
        for j= 2:sL;
            kk= zz+j-1;
            de(kk,:)= par(j,:);
        end
    end
end
sza= size(de,1);
dx1= Lx./(sza-1); 

for j= 1:sza
    de(j,1)= j-1;
end

%%%%%%%%%% Ouput %%%%%%%%%%%%%%%%%%%%%%%%%%%

save(vout, 'de', '-ascii');

%%%%%%%%%%  plots %%%%%%%%%%%%%%%%%%%%%%%%%%

szde= size(de,1);
szd2= size(de,2)
dnp= 7+mm;
ff = 1;
N= (szd2 - 1 - 2*ff)./dnp,      'Number of particle species'

nn(1)= 4+mm;            % density
vv(1)= 2;
vy(1)= 3+mm;
TT(1)= 5+mm;
qq(1)= 7+mm;
for ii=2:N
    nn(ii)= nn(ii-1) + dnp;
    vv(ii)= vv(ii-1) + dnp;
    vy(ii)= vy(ii-1) + dnp;
    TT(ii)= TT(ii-1) + dnp;
    qq(ii)= qq(ii-1) + dnp;
end
pt= szd2 - 1;
ee= szd2;

%%%%%%%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

colorset = [
	0.00  0.00  1.00
	0.00  0.50  0.00
	1.00  0.00  0.00
	0.00  0.75  0.75
	0.75  0.00  0.75
	0.75  0.75  0.00
	0.25  0.25  0.25
	0.75  0.25  0.25
	0.95  0.95  0.00
	0.25  0.25  0.75
	0.75  0.75  0.75
	0.00  1.00  0.00
	0.76  0.57  0.17
	0.54  0.63  0.22
	0.34  0.57  0.92
	1.00  0.10  0.60
	0.88  0.75  0.73
	0.10  0.49  0.47
	0.66  0.34  0.65
	0.99  0.41  0.23
];

%%%%%%% Densities %%%%%%%%%%%%%%%%%%%%%%%%%

fl=1;
figure(fl)
tmp=nn;

for ii=1:N
    plot(de(:,1),de(:,tmp(ii)), 'Color', colorset(ii,:))
    hold on;   
end

legend('1','2','3', '4', '5', '6', '7', '8', '9', '10','11','12','13','14','15');
xlabel('x [m]')
ylabel('n [m^-^3]')
axis tight

%%%%%%%%%%%% Parallel velocities %%%%%%%%%%%%%%

fl= fl + 1;
figure(fl)
tmp=vv;

for ii=1:N
    plot(de(:,1),de(:,tmp(ii)), 'Color', colorset(ii,:))
    hold on;   
end

legend('1','2','3', '4', '5', '6', '7', '8', '9', '10','11','12','13','14','15');

xlabel('x [m]')
ylabel('V_|_| [m/s]')
axis tight

%%%%%%%%%%%% Vy %%%%%%%%%%%%%%

fl= fl + 1;
figure(fl)
tmp=vy;

for ii=1:N
    plot(de(:,1),de(:,tmp(ii)), 'Color', colorset(ii,:))
    hold on;   
end

legend('1','2','3', '4', '5', '6', '7', '8', '9', '10','11','12','13','14','15');

xlabel('x [m]')
ylabel('V_y [m/s]')
axis tight

%%%%%%%%%%% Temperature  %%%%%%%%%%%%%%%

fl= fl + 1;
figure(fl)
tmp=TT;

for ii=1:N
    plot(de(:,1),de(:,tmp(ii)), 'Color', colorset(ii,:))
    hold on;   
end

legend('1','2','3', '4', '5', '6', '7', '8', '9', '10','11','12','13','14','15');

xlabel('x [m]')
ylabel('T [eV]')
axis tight

%%%%%%%%%%% Parallel heat fluxes %%%%%%%%%%%%%%%

fl= fl + 1;
figure(fl)
tmp=qq;

for ii=1:N
    plot(de(:,1),de(:,tmp(ii)), 'Color', colorset(ii,:))
    hold on;   
end

legend('1','2','3', '4', '5', '6', '7', '8', '9', '10','11','12','13','14','15');

xlabel('x [m]')
ylabel('q_x [W/m^2]')
axis tight

  
%%%%%%%%%%% Potential %%%%%%%%%%%%%%%%%%%%%%%%
    
    fl= fl + 1;
    figure(fl)
    
    plot(de(:,1),de(:,pt))
    xlabel('x [m]')
    ylabel('Pot. [V]')
    axis tight
    
    
%%%%%%%%%%% Electric field %%%%%%%%%%%%%%%%%%%%

    
    fl= fl + 1;
    figure(fl)
    
    plot(de(:,1),de(:,ee))
    xlabel('x [m]')
    ylabel('E_x [V/m]')
    axis tight

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

