clear
%script per la generazione di una processo di interarrivo
%con coefficiente di variazione compreso tra 1/2 e 3 (escluso)
tstart=clock;
R=100;   %realizzazioni
N=100;   %arrivi per realizzazione
tmi=1.2; %tempo medio di interarrivo
cv=2.5;  %cv=1 esponenziale, 1/2<=cv<1 erlang, 1<cv<3 iperesp
if cv<1/2
    disp('error')
    return;
else
    if cv<1
        type=1;
    else
        if cv==1
            type=0;
        else
            if cv<3
                type=2;
            else
                disp('error')
                return
            end
        end
    end
end
        
switch type
    case 0
        u=rand(R,N);
        ti=-tmi*log(1-u); %tempi di interarrivo
    case 1
        m(1,2)=tmi/2*(1+(2*cv-1)^.5);
        m(1,1)=tmi-m(1,2);
        for i=1:R
            u=rand(2,N);
            ti(i,:)=-m*log(1-u);
        end
    case 2
        l2=(tmi*(1+((cv-1)/2)^.5))^-1;
        l1=(2*tmi-1/l2)^-1;
        if l1<0
            l2=(tmi*(1-((cv-1)/2)^.5))^-1;
            l1=(2*tmi-1/l2)^-1;
        end
        if l1<0 || l2<0
            disp('error')
            return
        end
        lmax=max(l1,l2);
        lmin=min(l1,l2);
        xt=0:0.01:2*log(10)/lmax;
        ti=[];
        u=rand(R,N);
        for j=1:R
            i0=0;
            for i=1:N 
                [tia,tib]=find(exp(-xt*l1)+exp(-xt*l2)-2*(1-u(j,i))<0); 
                if tib 
                    i0=i0+1; 
                    ti(j,i0)=xt(min(tib)); 
                else
                    i0=i0+1;
                    ti(j,i0)=-log(2*(1-u(j,i)))/lmin;
                end
            end
        end
end
%keyboard
dt=0.01;
ta=cumsum(ti,2); %tempi di arrivo
tmin=ceil(max(min(ta'))); %ultimo primo tempo di arrivo
tmax=floor(min(max(ta'))); %primo ultimo tempo di arrivo
t=tmin+dt:dt:tmax-dt; %intervallo in cui ci sono arrivi in tutte le realizzazioni
for i=1:R
    if rem(i,50)==0, i, end
    for j=1:length(t) 
        tp=t(j)-ta(i,:);
        tla(i,j)=min(tp(tp>0));
        tna(i,j)=-max(tp(tp<0));
    end     
end
mtla=mean(tla); %tempo medio dall'ultimo arrivo
mtna=mean(tna); %tempo medio dal prossimo arrivo 
for i=1:3
    subplot(3,2,2*i)
    plot(t(t<3*tmin),tna(i,t<3*tmin),'k'); 
    set(gca, 'Xlim', [tmin,3*tmin]);
    set(gca, 'Ylim', [0,max(ceil(max(max(tna(1:3,t<3*tmin)))),ceil(max(max(tla(1:3,t<3*tmin)))))]);
    if i==1
        title('tempi per il prossimo arrivo (vita residua)');
    end
end

for i=1:3
    subplot(3,2,2*i-1)
    plot(t(t<3*tmin),tla(i,t<3*tmin),'k'); 
    set(gca, 'Xlim', [tmin,3*tmin]); 
    set(gca, 'Ylim', [0,max(ceil(max(max(tna(1:3,t<3*tmin)))),ceil(max(max(tla(1:3,t<3*tmin)))))]);
    if i==1
        title('tempi per dall"ultimo arrivo (vita)');
    end    
end
etime(clock,tstart)