/* Simulation der Epidemie Schistosomiasis */ /* Mathematische Spezifikation: Modell SNM, d.h. Barbour-Kafetzaki mit zusaetzlich moeglichen Todesfaellen und sofortige gekoppelte Geburt, konstante Menschenpopulation, F1 Poisson Ausdruck ist die Verteilung der Parasiten pro Mensch in gewissem Alter: 5/10/15/20/30/40/60/100. Das Programm ist folgendermassen aufgebaut: Programmkopf mit Importen Deklarationen Variablen- und Parameterbeschreibungen Initialisierungen und Voreinstellungen Hauptprogramm mit Zeitberechnung Ereignistypberechnung Neueinstellungen Programmende mit Ausdrucksberechnung Ausdruck */ /* --------------------------------------------------------------------------------------------------------- */ /* ------------------------------------------Programmkopf--------------------------------------------------- */ /* --------------------------------------------------------------------------------------------------------- */ /* --------------------------------------------Importe------------------------------------------------------ */ #include #include #include main() { /* int random(); bereits in stdlib deklariert, Funktion, die ich brauche fuer Uniform (0,1), Zahl zwischen 0 und 2E31-1 */ /* -----------------------------------------Deklarationen--------------------------------------------------- */ double abz,bp,dt,gas,glz,h,lambda,maxlife,maxsl,mlm,my,kappa,r1,r2,r3,r4,r5,r6; double hsts,sts,summe1,summe2,t,t1,theta; double prev[5001],prevn[5001],prevz[5001], d[101], d00[101], d02[101], d04[101], d06[101], d10[101], d15[101], d30[101], nschisto[101],thuman[200], ninf; int et,gam,i,j,k,kadenz,msm,no,shuman[200]; /* ------------------------------Variablen- und Parameterbeschreibungen------------------------------------- */ /* abz ist die Abbruchzeit bp ist der Belastungsparameter fuer die Aufteilung d[i] is the proportion of individuals with i worms over all age classes (distribution) dxy[i] is the proportion of individuals with i worms at age [xy, xy+1) dt ist die Zeit bis zum naechsten Ereignis et ist der Zustandsparameter bei der Frage nach dem Typ Ereignis gam ist die Gesamtanzahl Menschen (versehen mit Nummern Null bis (gam-1)) gas ist die momentane Gesamtanzahl Schistosomen im Programm (nur wenn gas gebraucht wird, nicht global) glz ist die aktuelle Gesamtlaufzeit h ist eine lokale Hilfsgroesse hsts Hilfe bei der sukzessiven Testsummation i ist ein nur lokal verwendeter Laufindex j ist ein nur lokal verwendeter Laufindex k ist ein nur lokal verwendeter Laufindex lambda ist die Kontaktrate der Menschen maxlife ist das maximal jemals erziehlte Alter eines Menschen (lebt eventuell nicht mehr) maxsl ist die maximale "Schistosomenload" (SL), welche ein Mensch im ganzen Programmablauf je hatte mlm ist die maximale Lebensmessung, welche noch registriert wird msm ist die maximale Schistosomenmessung, welche noch registriert wird (wegen Fakultaetsberechnungen (!)) my ist die Absterberate der Schistosomen im Menschen ninf is the sum for to calculate the overall prevalence (NonINFected) no ist das neue Opfer, der Mensch mit Nummer no wurde neu infiziert nschisto[p] ist die Anzahl Menschen mit p Schistosomen kappa ist die Absterberate der Menschen prev[i] is the proportion of infected at age [i, i+1) prevn[i] is the denominator for the calculations of prev[i] prevz[i] is the numerator for the calculations of prev[i] r1 bis r6 sind Zufallszahlen, welche alle mit Hilfe von random() gebildet werden, genaue Angaben bei Gebrauch shuman[t] ist die Anzahl Schistosomen in Person t summe1 ist die Gesamtrate mit der Schistosomen sterben summe2 ist die Gesamtrate mit der Menschen infiziert werden sts sukzessive Testsummation t ist eine lokale Hilfsgroesse t1 ist eine lokale Hilfsgroesse theta ist die durchschnittl. Anz. Parasiten, welche von einem Parasiten uebertragen werden (erfolgr. Infekt.) thuman[j] ist das Alter der Person j */ /* ----------------------------Initialisierungen und Voreinstellungen--------------------------------------- */ /* Die folgenden Initialisierungen und Voreinstellungen muessen konsistent und sinnvoll sein; es sind keine Tests eingebaut, d.h. gigo. */ /* Angaben ueber die Menschenpopulation: */ gam = 200; gas = 224; maxlife = 0; /* wird spaeter sofort korrigiert */ maxsl = 5; mlm = 5000; msm = 100; /* nicht vergroessern wegen Fakultaetsberechnungen (!) */ for (i=0; i<50; ++i) {shuman[i]=1;} for (i=50; i<60; ++i) {shuman[i]=2;} for (i=60; i<75; ++i) {shuman[i]=3;} for (i=75; i<91; ++i) {shuman[i]=4;} for (i=91; i<100; ++i) {shuman[i]=5;} for (i=100; imaxlife) { maxlife=thuman[i]; } } /* Zwischenmessung fuer die Praevalenz, small inacurracy: must be that dt << 1 */ ninf=ninf+nschisto[0]*dt; for(i=0;i<=maxlife;++i) { for(j=0;j=1) { prevz[i]=prevz[i]+dt; } } } } for(j=0;j<=gam;++j) { if ((thuman[j] >=0) && (thuman[j] < 1)) { d00[(shuman[j])]=(d00[(shuman[j])] + dt); } } for(j=0;j<=gam;++j) { if ((thuman[j] >=2) && (thuman[j] < 3)) { d02[(shuman[j])]=(d02[(shuman[j])] + dt); } } for(j=0;j<=gam;++j) { if ((thuman[j] >=4) && (thuman[j] < 5)) { d04[(shuman[j])]=(d04[(shuman[j])] + dt); } } for(j=0;j<=gam;++j) { if ((thuman[j] >=6) && (thuman[j] < 7)) { d06[(shuman[j])]=(d06[(shuman[j])] + dt); } } for(j=0;j<=gam;++j) { if ((thuman[j] >=10) && (thuman[j] < 11)) { d10[(shuman[j])]=(d10[(shuman[j])] + dt); } } for(j=0;j<=gam;++j) { if ((thuman[j] >=15) && (thuman[j] < 16)) { d15[(shuman[j])]=(d15[(shuman[j])] + dt); } } for(j=0;j<=gam;++j) { if ((thuman[j] >=30) && (thuman[j] < 31)) { d30[(shuman[j])]=(d30[(shuman[j])] + dt); } } for(i=0;i<=msm;++i) { d[i]=(d[i] + (nschisto[i]*dt)); } /* ---------------------------------------Ereignistypberechnung--------------------------------------------- */ r2 = ((((double)random())+1)/2147483650.0); if (r2< (summe1/bp)) et=0; /* Ein Schistosome stirbt */ else if ((r2 >= (summe1/bp)) && (r2 < ((summe1+summe2)/bp))) et=1; /* Eine Infektion passiert */ else et=2; /* Ein Mensch stirbt */ /* -----------------------------------------Neueinstellungen------------------------------------------------ */ /* Erste Veraenderungsmoeglichkeit: ein Schisto starb(et=0) */ if (et==0) { /* Progr.strat.: Schistos aufgereiht in Menschenreihenfolge */ gas=0; for (i=1; i<=maxsl; ++i) { gas=(gas+(i*nschisto[i])); } r3=(((((double)random())+1)/2147483650.0)*gas); sts=0; for (i=0; i< gam; ++i) { if((sts <= r3) && (r3 < (sts+shuman[i]))) { /* Sch. aus Mensch i starb */ nschisto[shuman[i]]=(nschisto[shuman[i]]-1);/* ein dekrement */ nschisto[(shuman[i]-1)]=(nschisto[(shuman[i]-1)]+1); /* incr. */ sts=(sts+shuman[i]); /* muss hier erhoehen falls nur 1 Schisto */ shuman[i]=(shuman[i]-1); /* Person i hat einen Schisto weniger */ } sts=(sts+shuman[i]); } } /* Zweite Veraenderungsmoeglichkeit: eine Infektion passierte(et=1), ev. Nullinfektion wegen bp-Berechnung */ else if (et==1) { /* Welcher Mensch wird infiziert? */ sts=0; r4 = (((((double)random())+1)/2147483650.0)*nschisto[0]); no=0; while(sts<=r4) { if(shuman[no]==0) { ++sts; } ++no; } no=(no-1); /* Mensch no wurde infiziert */ /* Mit wievielen Schistosomen wird Mensch no infiziert? */ sts=0; for (i=0;i<=msm;++i) { hsts=0; if(i==0) { t1=1; } else { t1=t1*i; /* t1 ist Fakultaet von i */ } for (j=1;j<=maxsl;++j) { h=j*theta; t=i; hsts=hsts+nschisto[j]*(exp((-h)))*pow(h,t)/t1; } sts=sts+hsts; } /* sts ist nun proportional zum Gesamtmass der Infektionsanteile */ r5 = (((((double)random())+1)/2147483650.0)*sts); sts=0; i=0; while(r5>sts) { hsts=0; if(i==0) { t1=1; } else { t1=t1*i; /* t1 ist Fakultaet von i */ } for (j=1;j<=maxsl;++j) { h=j*theta; t=i; hsts=hsts+nschisto[j]*(exp((-h)))*pow(h,t)/t1; } sts=sts+hsts; ++i; } k=i-1; /* Person no wird mit k Schistosomen infiziert */ nschisto[k]=(nschisto[k]+1); nschisto[0]=(nschisto[0]-1); /* ev. Nullinfektion */ shuman[no]=k; if (k>maxsl) { maxsl=k; /* Die max. SchistoLoad wird erhoeht */ } } /* Schluss der else if (et==1) */ /* Dritte Veraenderungsmoeglichkeit: ein Mensch stirbt(et=2) */ else if (et==2) { /* Progr.strat.: Menschen aufgereiht */ sts=0; r6 = (((((double)random())+1)/2147483650.0)*gam); for(i=0;i