function [Type,M_Evol,IFS_Evol,P_macro,Event_macro_Evol,firing_transition,timer_macro_event,tau,Q_Evol,nc_nd_qc_qd]=simulator_HPN(Pre,Post,M,vel,v,D,s,alpha,time_stop,simulation_type,it)
%
% simulator_HPN : This function simulates a Hybrid Petri Net
% *********************************************************************
% ## SYNTAX ##
%
% [Type,M_Evol,IFS_Evol,P_macro,Event_macro_Evol,firing_transition,
% timer_macro_event,tau,Q_Evol,nc_nd_qc_qd]=simulator_HPN(Pre,Post,
% M,vel,v,D,s,alpha,time_stop,simulation_type,it)
%
% Input arguments:
%
% Pre: Pre-incidence matrix
% |Pre_CC NaN Pre_CD |
% |NaN NaN NaN NaN NaN|
% |Pre_DC NaN Pre_DD |
%
% where Pre_XY --> X=C: continuous places
% X=D: discrete places
% Y=C: continuous transitions
% Y=D: discrete transitions
%
% Post: Post-incidence matrix
% |Post_CC NaN Post_CD|
% |NaN NaN NaN NaN NaN|
% |Post_DC NaN Post_DD|
%
% where Post_XY --> X=C: continuous places
% X=D: discrete places
% Y=C: continuous transitions
% Y=D: discrete transitions
%
% M: Initial marking
% [M0_C;
% M0_D]
%
% vel: matrix (qc x 2) that specifies, for each continuous transition,
% the mfs and MFS
%
% |msf(1) MFS(1) |
% | ... ... |
% |msf(qc) MFS(qc)|
%
% v: qd-dimensional row vector that specifies the timing structure
% of discrete transitions
%
% [v(1),...,v(qd)]^(T)
%
% Deterministic: 1, Exp: 2, Unif: 3, Poiss: 4, Rayl: 5,Wbl:6 ,
% Beta: 7, Chi2: 8, Ev: 9, F: 10, gam: 11, gev: 12, gp: 13,
% logn: 14, ncf: 15, nct: 16, ncx2:17, norm: 18, t: 19.
% Consult "help random" to get more information
%
% D: (qd x 3) matrix that associates to each discrete transition
% a row vector of length 3. If the transition is deterministic
% the first element of the row is equal to the time delay of
% transition; if the transition is stocastic, the elements of
% the row specifies the parameters of the corresponding
% distribution function (up to 3)
%
% |par_1(1) par_2(1) par_3(1) |
% | ..... ..... ..... |
% |par_1(qd) par_2(qd) par_3(qd)|
%
% s: qd-dimensional row vector that keeps track of the number of
% servers associated to discrete transitions
%
% [s(1) ... s(qd)]
%
% 0: infinite-server semantics, k: k-server semantics
%
% alpha: specifies the conflict resolution policy among discrete
% transitions
%
% |alpha_1(1) ..... alpha_1(qd)| <--transitions weight
% |alpha_2(1) ..... alpha_2(qd)| <--Priority level
%
% In the simulator several conflict resolution policies can be
% choosen:
% - Fire according to the weight of each transition
% - Fire according to the index of each transition
% - Fire according to the priority level and to the weight of
% each transition
%
% time_stop: is equal to the time length of simulation
%
% simulation_type: specifies the simulation mode
% - 2: Generates on screen the evolution graph (the simulation is
% carried out step-by-step)
% - 1: Generates on screen the evolution graph (the simulation is
% carried out without interruptions until time_stop is
% reached)
% - 0: No intermediate result is shown but only array Evol
% is created
%
% it: language parameter
% - 1: Italian
% - 0: English
%
% When function simulator_HPN is executed a qc-dimensional row
% vector J is required from the workspace, that associates to each
% continuous transition a weight
%
% J=[J(1) ... J(qc)]
%
% J*V is the linear cost function that should be maximized at each
% MP to compute the IFS vector V
%
% *********************************************************************
% Output arguments:
%
% (In the following k is the number of ME that occur during the
% simulation run)
%
% Type: 1(2,3) if the net is continuous(discrete,hybrid)
%
% M_Evol: matrix (n x (k+1)) that keeps track of the marking of the
% net during all the evolution. An n-dimensional column
% vector is associated to the initial time instant and to
% all the time instants in which a different ME occurs,
% each one representing the corresponding value of the
% marking at that time instant
%
% IFS_Evol: matrix (qc x (k+1)) that keeps track of the IFS vector
% during all the evolution. A qc-dimensional column vector
% is associated to the initial configuration and to the
% end of each ME
%
% P_macro: (k+1)-dimensional row vector that keeps track of the ME
% caused by continuous places
% The r-th entry may take the following values:
% * i : the (r+1)-th ME is caused by continuous place pi
% * NaN : the (r+1)-th ME is caused by a discrete transition
%
% Event_macro_Evol: (k+1)-dimensional row vector that keeps track
% of the ME caused by continuous places
% The r-th entry may take the following values:
% * 0 : the continuous place gets empty
% * +1 : the continuous place enables a new discrete
% transition
% * -1 : the continuous place disables a discrete
% transition
% * NaN : the r-th ME is due to the firing of a
% discrete transition
%
% firing_transition: (k+1)-dimensional row vector that keeps track
% of the discrete transitions that have fired
% during the evolution
% The r-th entry may take the following values:
% * i : the (r+1)-th ME is caused by discrete transition ti
% * 0 : the (r+1)-th ME is caused by a continuous place
%
% timer_macro_event: (k+1)-dimensional row vector that keeps into
% memory the length of ME.
%
% tau: total time of simulation
%
% Q_Evol: array of cell ((k+1) x nd) whose generic (r+1,j)-th
% entry specifies the clocks of transition tj at the end
% of the r-th MP
%
% nc_nd_qc_qd: 4-dimensional row vector equal to [nc nd qc qd]
%
% *********************************************************************
% Other sintax
%
% 1) [Type,M_Evol,IFS_Evol,P_macro,Event_macro_Evol,firing_transition,
% timer_macro_event,tau,Q_Evol,nc_nd_qc_qd]=simulator_HPN(net)
%
% net={Pre,Post,M0,vel,v,D,s,alpha}
% ==>
