Função Calcula o valor de d (deslocamento no tempo) para qual o fator HIC (Head Injury Criterion) é máximo

// Calcula o valor de d (deslocamento no tempo) para o qual o fator HIC 

// (Head Injury Criterion) é máximo, usando as fórmulas para aceleração

// sem airbag (a_s) e com airbag (a_c), a partir de tempos iniciais t_0_s 

// e t_0_c (em milissegundos)

a_s='16400/((x-68)**2+400)+1480/((x-93)**2+18)'

a_c='22000/((x-74)**2+500)'

t_0_s=63

t_0_c=56

d=(3:1:72)'

max_hic_s=0

max_hic_c=0

clear hic_c hic_s

for i=1:max(size(d))

    hic_s(i)=1e-3*d(i)*(1/d(i)*simpson_exat_crescente(t_0_s,t_0_s+d(i),8,a_s))**2.5

    hic_c(i)=1e-3*d(i)*(1/d(i)*simpson_exat_crescente(t_0_c,t_0_c+d(i),8,a_c))**2.5

    if hic_s(i)>max_hic_s then

        max_hic_s = hic_s(i)

        d_s = d(i)

    end

    if hic_c(i)>max_hic_c then

        max_hic_c = hic_c(i)

        d_c = d(i)

    end    

end

mprintf('sem airbag: valor máximo HIC=%g em d=%g\n',max_hic_s,d_s)

mprintf('com airbag: valor máximo HIC=%g em d=%g\n',max_hic_c,d_c)

clf

plot(d,hic_c,'r-',d,hic_s,'g--')

set(gca(),"grid",[1 1])

xtitle('d x HIC')

xlabel('d')

ylabel('HIC')