real S,I,R,C,D real gama,mu,beta real t,dt,tm,pm real ni,s0,i0,r0,tt real di,mu1,t1,III integer ii1,j00,nmm integer j1 parameter(nmm=1000) real tempos(nmm),dados(nmm) character*100 file100,file200 real s10(0:nmm),r10(0:nmm),i10(0:nmm) real s20(nmm,nmm),r20(nmm,nmm),i20(nmm,nmm),minn real diff(nmm,nmm),betaa(nmm,nmm) real beta100(nmm),gamaa(nmm,nmm),gama100(nmm) integer nd,ii,nbb,ntt,itt,i00,idd real dadosm(0:nmm),mort(0:nmm),mmm100 real aux1000,r1,r2,r3 real R01(0:nmm),R02(0:nmm),R03(0:nmm),R04(0:nmm) integer i000,j000 integer ii000,jj000 real beta200(nmm),gama200(nmm) real beta300(nmm),gama300(nmm) real beta400(nmm),gama400(nmm) real RRR(0:nmm),errom i000=1 j000=1 ii000=1 jj000=1 dt=0.00001 tm=71. pm=220000000. ! ENTRA COM O TAMANHO DOS DADOS DE ESTRADA (o numero de dias depois de 22 de março) write(*,*)"n" read(*,*)nd ! Ler as tabelas de casos confirmados e casos letais 1933 format(a100) ! write(*,*)"dados de entrada casos confirmados" ! read(*,1933)file100 file100="tabelabrasil.txt2" open(12,file=file100) ! write(*,*)"dados de entrada casos letais" !read(*,1933)file200 file200="mortesbrasil.txt2" open(122,file=file200) do ii=1,nd read(12,*)tempos(ii),dados(ii) read(122,*)tempos(ii),dadosm(ii) enddo idd=1 di=tempos(idd) ni=dados(idd) ! Inicia S, I e R para t=0 s10(0)=1-2.*ni/pm r10(0)=ni/pm i10(0)=ni/pm !Criar nomes para os arquivos de saida open(15,file="BAinfectados.dat") ! Arquivo principal com os resultados de I ao longo do tempo open(110,file="BAR0medioxtempo.dat") open(1100,file="BAR0xtempo.dat") open(17,file="BAdadosreaiscasosconfirmados.dat") open(177,file="BAdadosreaiscasosletais.dat") open(220,file="BAcalosletais.dat") !gama=0.09 mu=0.015 mu1=0.035 ntt=0 t1=0. itt=0 ! Inicia a simulacao partindo de t=0 (consideramos 22 de março no Brasil ) até o dia final (o dia presente) do temp=1.0e0,float(nd-1),1.0e0 idd=idd+1 ! Roda gama e beta dentro da faixa escolhida naa=0 do gama=0.08,0.115,0.005 naa=naa+1 nbb=0 do beta=0.08,0.5,0.005 ! beta nbb=nbb+1 s0=s10(itt) r0=r10(itt) i0=i10(itt) do t=dt+t1,t1+1.,dt !tempo !if (t.gt.50) beta=0.2 ! if (t.gt.75) beta=0.3 !Integrar o modelo SIR S=s0-dt*beta*i0*s0 I=i0+dt*(beta*i0*s0-gama*i0) R=r0+dt*gama*i0 s0=S i0=I r0=R enddo !tempo s20(nbb,naa)=S r20(nbb,naa)=R i20(nbb,naa)=I III=I*pm betaa(nbb,naa)=beta gamaa(nbb,naa)=gama diff(nbb,naa)=abs(III-dados(idd)) enddo ! beta enddo ! gamma minn=1e+20 do ii1=1,naa do ii=1,nbb if (diff(ii,ii1).lt.minn)then minn=diff(ii,ii1) ! Escolhe qual valores de gama e beta ! fornecem valores de I proximo dos dados reais !Escolhemos 3 conjuntos de valores de gama,beta ! aqueles 3 que ficam mais proximos ! o e é medido atraves d euma média destes 3 valores iii000=ii000 jjj000=jj000 ii000=i000 jj000=j000 i000=i00 j000=j00 i00=ii j00=ii1 endif enddo enddo itt=itt+1 s10(itt)=s20(i00,j00) r10(itt)=r20(i00,j00) i10(itt)=i20(i00,j00) beta100(itt)=betaa(i00,j00) gama100(itt)=gamaa(i00,j00) beta200(itt)=betaa(i000,j000) gama200(itt)=gamaa(i000,j000) beta300(itt)=betaa(ii000,jj000) gama300(itt)=gamaa(ii000,jj000) beta400(itt)=betaa(iii000,jjj000) gama400(itt)=gamaa(iii000,jjj000) !write(*,*)itt write(15,*)temp,i10(itt)*pm !write(,dados(idd) write(17,*)temp,dados(idd) write(177,*)temp,dadosm(idd) !write(110,*)temp,beta100(itt)/gama100(itt),beta100(itt),gama100(itt) mort(itt)=dadosm(idd)/dados(idd) write(220,*)temp,mort(itt)*i10(itt)*pm t1=t1+1.0e0 R01(itt)=beta100(itt)/gama100(itt) R02(itt)=beta200(itt)/gama200(itt) R03(itt)=beta300(itt)/gama300(itt) R04(itt)=beta400(itt)/gama400(itt) enddo ! tempo !R0 ao longo do tempo do ii=1,itt RRR(ii)=R01(ii)+R02(ii)+R03(ii)+R04(ii) RRR(ii)=RRR(ii)*0.25 enddo errom=0. do ii=1,itt write(1100,*)ii,RRR(ii),beta100(itt),gama100(itt) errom=errom+(RRR(ii)-R01(ii))**2. errom=errom+(RRR(ii)-R02(ii))**2. errom=errom+(RRR(ii)-R03(ii))**2. errom=errom+(RRR(ii)-R04(ii))**2. enddo errom=errom/(float(4*itt)) write(*,*)errom/RRR(itt) do ii=4,itt-3,4 r1=RRR(ii)+RRR(ii-1)+RRR(ii-2)+RRR(ii-3) r1=r1*0.25 r2=beta100(ii)+beta100(ii-1) r2=r2+beta100(ii-2)+beta100(ii-3) r2=r2*0.25 r3=gama100(ii)+gama100(ii-1) r3=r3+gama100(ii-2)+gama100(ii-3) r3=r3*0.25 write(110,*)ii,r1,r2,r3 enddo ii=itt r1=RRR(ii)+RRR(ii-1)+RRR(ii-2)+RRR(ii-3) r1=r1*0.25 r2=beta100(ii)+beta100(ii-1) r2=r2+beta100(ii-2)+beta100(ii-3) r2=r2*0.25 r3=gama100(ii)+gama100(ii-1) r3=r3+gama100(ii-2)+gama100(ii-3) r3=r3*0.25 write(110,*)ii,r1,r2,r3 write(*,*)r1,r2,r3 gama=r3 beta=0.999*r3*r1 mmm100=(mort(itt)+mort(itt-1)+mort(itt-2)) mmm100=mmm100+mort(itt-3)+mort(itt-4) mmm100=mmm100/5. ! Aqui começa a previsão para 10 dias ! Usando o ultimo valor de beta e gama ! calculado como sendo a média da ultima semana de dados t1=float(nd-1) do temp=float(nd),float(nd)+10.,1.0e0 !tempo 2 s0=s10(itt) r0=r10(itt) i0=i10(itt) do t=dt+t1,t1+1.,dt !tempo X !if (t.gt.50) beta=0.2 ! if (t.gt.75) beta=0.3 !integra as equações S=s0-dt*beta*i0*s0 I=i0+dt*(beta*i0*s0-gama*i0) R=r0+dt*gama*i0 s0=S i0=I r0=R enddo !tempo X itt=itt+1 s10(itt)=s0 r10(itt)=r0 i10(itt)=i0 t1=t1+1.0e0 write(15,*)temp,i10(itt)*pm write(220,*)temp,mmm100*i10(itt)*pm write(110,*)temp,beta/gama !if (temp.gt.100) then !if (beta.gt.0.09) beta=0.991*beta !endif !O beta o Brasil tem uma leva diminuição ao longo das semanas ! mantivemos esta diminuição usando este truque de multiplicar diariamente !o beta pro uma pequena constante ligeiramente menor que 1 beta=0.999*beta mmm100=0.9999*mmm100 enddo ! loop do temp 2 end