%% this program analyzes data from the "/data/" folder
clc;clear all;
nstep=500;
Re=6371e3;h=1e-3;np=5000;
vp=[];x=[];y=[];z=[];vper=[];

for i=1:nstep
    t(i)=i*h;
    
fname=['vp_' num2str(i) '.dat'];
a=load(fname);
vp=[vp;a];
lk=[];
lk=find(isnan(a(:)));
loss_count(i)=length(lk);

fname=['vper_' num2str(i) '.dat'];
a=load(fname);
vper=[vper;a];

fname=['x_' num2str(i) '.dat'];
a=load(fname);
x=[x;a];

fname=['y_' num2str(i) '.dat'];
a=load(fname);
y=[y;a];

fname=['z_' num2str(i) '.dat'];
a=load(fname);
z=[z;a];

end
[phi,lam,r] = cart2sph(x,y,z);
phi=phi.*180/pi;
lam=lam.*180/pi;

figure(1)
hold on
plot((r-Re)./1e3,lam,'.')
xlabel('height [km]');
ylabel('\lambda [deg]');

rr=(r-Re)./1e3;
[l1 l2]=size(vp);
vp=reshape(vp,[l1*l2,1]);
r=reshape(r,[l1*l2,1]);
lam=reshape(lam,[l1*l2,1]);
phi=reshape(phi,[l1*l2,1]);
vper=reshape(vper,[l1*l2,1]);
lk=find(isnan(vp(:)));
vp(lk)=[];vper(lk)=[];
r(lk)=[];lam(lk)=[];phi(lk)=[];

lk=length(vp);
l1=0;
for j=1:lk
    if abs(lam(j))<6
        l1=l1+1;
        lam1(l1)=lam(j);
        r1(l1)=r(j);
        vp1(l1)=vp(j);
        phi1(l1)=phi(j);
        vper1(l1)=vper(j);
        al_cal(l1)=acosd(vp(j)/sqrt(vp(j)*vp(j)+vper(j)*vper(j)));
    end
end
%% plotting

figure(2)
subplot(2,2,1)
histogram(vp1./1e3)
xlabel('v_{par} [km/s]')
ylabel('Occurrence')

subplot(2,2,2)
histogram(vper1./1e3)
xlabel('v_{per} [km/s]')
ylabel('Occurrence')

subplot(2,2,4);
histogram(al_cal);
xlabel('\alpha [deg]')
ylabel('Occurrence')

subplot(2,2,3);
plot(t,loss_count.*100/np,'.-');
xlabel('time [seconds]')
ylabel('Particle loss [%]')