clear;
saveDirectory='C:\Users\bacaer\Desktop\Nicolas\2010HIVTBarticles\Revise\';
saveDirectory2='C:\Users\bacaer\Desktop\Nicolas\2010Londres\Expose\';
loadDirectory='C:\Users\bacaer\Desktop\';
universal=0; //0 means no, 1 means yes
showassumptions=0; //same
savefig=0; //same
show='HIV';
//show='DEATH';
//show='PYRAMID';
//************** assumptions **********************
xmax=100;
ymax=xmax;
t0=1986;
liste=(20/1152)*[1152 1150 1147 1142 1136 1124 1109 1093 1074 1053 1031 1031 1030 1028 1025 1022	1018 1013 1008 1004 1000 993 987 978 963 955 943 930 916 906 902 898 893 887 880 874 868 861 854  846 840 834 828 821 814 808 803 797 791 785 780 774 769 764 759]; //births between 1986 and 2020 (US Census Bureau)
function z=B(t)
  z=liste(t-t0+1)
endfunction
tmax=2020;
Sm=zeros(xmax,1); Im=zeros(xmax,ymax); Rm=zeros(xmax,ymax); SSm=zeros(xmax,1); IIm=zeros(xmax,ymax); RRm=zeros(xmax,ymax);  //males
Sf=zeros(xmax,1); If=zeros(xmax,ymax); Rf=zeros(xmax,ymax); SSf=zeros(xmax,1); IIf=zeros(xmax,ymax); RRf=zeros(xmax,ymax); //females
df=0.001*linear_interpn([1:xmax],[1 2 5 10:5:100],[40 2 1 0.5 0.6 2.0 2.5 3.1 3.9 4.5 5.0 7.5 10.0 15 22 29 36 43 50 57 64 71])';  //female death rate
dm=0.001*linear_interpn([1:xmax],[1 2 5 10:5:100],[50 2 1 0.6 0.7 2.3 3.5 4.2 5.0 6.5 8.5 12.0 17.0 26 37 48 59 70 81 92 103 114])';  //male death rate
debut=14;  //age of female sexual debut
choix=zeros(xmax,xmax); ym=zeros(1,xmax); yave=zeros(1,xmax); yM=zeros(1,xmax);
//choice of male sexual partner
for x=debut:xmax,
  ym(x)=linear_interpn(x,[debut 40 100],[debut 22 84]); //minimum
  yave(x)=linear_interpn(x,[debut 30 100],[debut+1 30 95]); //modal
  yM(x)=linear_interpn(x,[debut 30 100],[debut+15 65 100]); //maximum
  choix(x,ym(x):yM(x))=linear_interpn(ym(x):yM(x),[ym(x) yave(x) yM(x)],[0 1 0]);
  choix(x,ym(x):yM(x))=choix(x,ym(x):yM(x))/sum(choix(x,ym(x):yM(x))); //probability distribution
end;
underreport=2; //underreporting of male sexual partners
rate=underreport*linear_interpn([1:xmax],[0 debut 15 25 35 45 60 100],[0 0 0.6 0.3 0.15 0.05 0.025 0.02])'; //rate of change of sexual partner
pm=0.5; //transmission probability female to male per partnership
pf=0.85;  //male to female 
overreport=1;  //overreporting of condom use
epsi=0.01; //relative infectivity on ART
AIDS=zeros(xmax,ymax);
//survivalAIDS=exp(-0.00494*([0:ymax].^2.25)); //median 9 years
survivalAIDS=exp(-0.0039*([0:ymax].^2.25)); //median 10 years
//survivalAIDS=exp(-0.00315*([0:ymax].^2.25)); //median 11 years
AIDSdeath=(survivalAIDS(1:ymax)-survivalAIDS(2:ymax+1))./survivalAIDS(1:ymax);
for x=1:xmax,
  if x<4 then
    AIDS(x,:)=0.4*(1-x/4)+(x/4)*AIDSdeath; //extra mortality for children with AIDS
  else
    AIDS(x,:)=AIDSdeath;
  end;
end;
ARTdeath=AIDS/2; //ART death rate
//non-normalized age distribution of mothers age
weight=0.001*linear_interpn([1:xmax],[1 14.9 15 19.9 20 24.9 25 29.9 30 34.9 35 39.9 40 44.9 45 49.9 50 xmax],[0 0 70.7 70.7 139 139 141.8 141.8 105.6 105.6 67.4 67.4 27.1 27.1 8.8 8.8 0 0])';
function z=condom(t)
  maxicondom=linear_interpn(t,[t0 1987 1995 2007 tmax],[0 0 0.1 0.8 0.8])
  decline=linear_interpn(t,[t0 1998 2008 tmax],[0.03 0.03 0.01 0.01])
  z=zeros(xmax,1)
  z(debut:xmax)=overreport*maxicondom*exp(-decline*([debut:xmax]'-debut).^1.2)//condom use decline with age
endfunction
drop=0.015; //ART drop out rate
chance=0.4; //probability of ending up on ART under current program
ART_intervention1=AIDS*chance/(1-chance);
ART_intervention2=0.2;   //annual proportion starting ART if universal treatment 
function z=ART(t)
  ART_access1=linear_interpn(t,[t0 2003 2007 tmax],[0 0 1 1])
  ART_access2=universal*linear_interpn(t,[t0 2010 2012 tmax],[0 0 1 1])
  z=ART_access1*ART_intervention1+ART_access2*ART_intervention2
endfunction
function z=M2C(t) //transmission probability mother to child as a function of time
  z=linear_interpn([t],[t0 2003 2008 tmax],[0.4 0.4 0.4*0.25+0.05*0.75 0.4*0.25+0.05*0.75])
endfunction
oo=ones(5,1);
//male age pyramid in 1986 US Census bureau by 5 year age groups (only for showing)
Sm0=(B(t0)/2553/1.07)*[2553*oo; 2202*oo; 1987*oo; 1763*oo; 1590*oo; 1434*oo; 1219*oo; 993*oo; 812*oo; 678*oo; 547*oo; 432*oo; 324*oo; 231*oo; 154*oo; 91*oo; 43*oo; 14*oo; 2*oo; 0.3*oo]; 
//smoothed age pyramid used in the simulation
Sm=(B(t0)/2553/1.07)*linear_interpn([0:xmax-1],[0 1 6 12:5:102],[1.07*2553 2553 2202 1987 1763 1590 1434 1219 993 812 678 547 432 324 231 154 91 43 14 2 0.3 0])';
Sf0=(B(t0)/2545/1.07)*[2545*oo; 2198*oo; 1995*oo; 1802*oo; 1651*oo; 1517*oo; 1304*oo; 1089*oo; 883*oo; 755*oo; 607*oo; 505*oo; 383*oo; 291*oo; 217*oo; 147*oo; 82*oo; 33*oo; 8*oo; 1*oo];
Sf=(B(t0)/2545/1.07)*linear_interpn([0:xmax-1],[0 1 6 12:5:102],[1.07*2545 2545 2198 1995 1802 1651 1517 1304 1089 883 755 607 505 383 291 217 147 82 33 8 1 0])';
//******************** age pyramid data **************
function zz=showpyramid()
  zz=0; clf(); xset("thickness",4); xset("font size",6);
  for i=1:17,
    x=[1:5]+5*(i-1);
    Pyramidm(i)=sum(Pm(x))/5;
    Pyramidf(i)=sum(Pf(x))/5;
  end;
  rrect=[-25,0,25,80]; nnax=[0,3,0,5];
//  plot2d2([-Pyramidm Pyramidf],[2.5:5:82.5]'*ones(1,2),style=[1 1],rect=rrect,nax=nnax);
  plot2d2([-Pm Pf],[0:xmax-1]'*ones(1,2),style=[1 1],rect=rrect,nax=nnax);
  plot2d([0 0],[0 80],rect=rrect,nax=nnax);
  xstring(15,75,string(t)); xstring(-20,50,"MEN"); xstring(10,50,"WOMEN"); xstring(-24,75,"AGE");
    if t==1990 then
     plot2d(-Pyramidm(1)*[2699128 2422207 2141571 1936475 1704933 1537846 1375270 1150019 929826 756410 622263 487078 371303 262781 172697 101621 48007 15863 3223 357]/2699128,[2.5:5:97.5],style=-9,rect=rrect,nax=nnax);
    plot2d(Pyramidf(1)*[2684387 2423851 2141283 1952517 1755993 1613175 1469895 1245508 1032798 833193 710514 558366 457793 329331 236141 158240 88029 35737 9209 1316]/2684387,[2.5:5:97.5],style=-9,rect=rrect,nax=nnax);
    end;
     if t==1995 then
     plot2d(-Pyramidm(1)*[2614638 2660223 2412492 2121742 1894307 1659325 1496587 1335764 1108490 884750 707953 570631 432274 313352 203841 116756 54767 18092 3664 410]/2614638,[2.5:5:97.5],style=-9,rect=rrect,nax=nnax);
    plot2d(Pyramidf(1)*[2599837 2646945 2409747 2115493 1910701 1720716 1581988 1441888 1215605 998333 796138 670853 515860 406652 272524 171213 93949 38108 9782 1401]/2599837,[2.5:5:97.5],style=-9,rect=rrect,nax=nnax);
    end;
   if t==2000 then
     plot2d(-Pyramidm(1)*[2445814 2560634 2633243 2387921 2061216 1809314 1572366 1411341 1255354 1032921 813997 639575 501458 362951 244138 139869 64608 21527 4511 514]/2445814,[2.5:5:97.5],style=-9,rect=rrect,nax=nnax);
    plot2d(Pyramidf(1)*[2432565 2564711 2633417 2380628 2044320 1830479 1644042 1519941 1391941 1169661 952659 750344 620155 460052 340027 203528 105612 43334 11532 1713]/2432565,[2.5:5:97.5],style=-9,rect=rrect,nax=nnax);
    end;
    if t==2005 then
     plot2d(-Pyramidm(1)*[2400652 2393233 2544612 2764641 2535108 2077890 1709531 1423408 1262752 1117155 912480 710936 546839 413134 280353 168745 79728 26929 5824 709]/2400652,[2.5:5:97.5],style=-9,rect=rrect,nax=nnax);
    plot2d(Pyramidf(1)*[2381255 2387101 2551350 2683348 2396257 1917306 1606211 1456503 1387651 1294646 1092013 882700 683323 547148 383682 257544 131108 52403 14640 2357]/2381255,[2.5:5:97.5],style=-9,rect=rrect,nax=nnax);
    end;
endfunction
//******************** mortality data ***********************
function zz=showdeath()
  zz=0; clf(); xset("thickness",4); xset("font size",6); rrect=[-50,0,50,70]; nnax=[4,3,0,8];
  for i=1:14,
    x=[1:5]+5*(i-1);
    MaleMortality(i)=(sum(dm(x).*(Sm(x)+TIm(x)+TRm(x)))+sum(Im(x,:).*AIDS(x,:))+sum(Rm(x,:).*ARTdeath(x,:)))/sum(Pm(x));
    FemaleMortality(i)=(sum(df(x).*(Sf(x)+TIm(x)+TRf(x)))+sum(If(x,:).*AIDS(x,:))+sum(Rf(x,:).*ARTdeath(x,:)))/sum(Pf(x));    
  end;
  plot2d2(1000*[-[MaleMortality(1); MaleMortality] [FemaleMortality(1); FemaleMortality]],[0:5:70]'*ones(1,2),style=[1 1],rect=rrect,nax=nnax);
  plot2d([0 0],[0 70],rect=rrect,nax=nnax);
  xstring(35,65,string(t)); xstring(-30,15,"MEN"); xstring(20,15,"WOMEN"); xstring(31,2,"/1000"); xstring(-15,65,"AGE");
  if t==1997 then // deaths from report p169, population from IDB
    plot2d(-1000*[45730 4030 3320 5180 7530 8580 10160 11720 12280 13000 14110 15350 16420]./[2519268 2641663 2528371 2204202 1946637 1711538 1523157 1374450 1174028 940544 744489 596216 458958],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);  
    plot2d(1000*[37560 2440 1850 3540 7120 7960 7160 7040 6910 7400 8530 10120 11710]./[2515141 2638488 2531958 2201467 1956715 1768463 1612416 1487179 1294070 1065345 850952 700535 555735],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(-1000*[17025 1704 1546 3775 8166 10908 11808 11942 11742 12200 11270 12616 11165]./[2519268 2641663 2528371 2204202 1946637 1711538 1523157 1374450 1174028 940544 744489 596216 458958],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);  
    plot2d(1000*[15189 1252 1188 2470 5403 7359 7134 6813 6365 6334 6225 7914 9283]./[2515141 2638488 2531958 2201467 1956715 1768463 1612416 1487179 1294070 1065345 850952 700535 555735],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
  end;
  if t==1998 then
    plot2d(-1000*[47690 4100 3220 5310 8350 10830 12960 14580 14780 14620 15200 15860 16860]./[2486301 2624041 2576391 2257640 1979562 1745300 1539538 1390842 1207480 972982 767491 610495 474311],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);  
    plot2d(1000*[39960 2570 1810 3870 9090 11040 9240 8550 7930 7970 8900 10340 12000]./[2481005 2622962 2579373 2254785 1982414 1794017 1626601 1502957 1333402 1100403 885109 714669 579468],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(-1000*[19780 1778 1691 4101 8781 13047 14339 14552 13887 14143 12962 13884 12399]./[2486301 2624041 2576391 2257640 1979562 1745300 1539538 1390842 1207480 972982 767491 610495 474311],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);  
    plot2d(1000*[17733 1435 1284 2895 6864 9755 9641 8854 7884 7637 7189 8860 9985]./[2481005 2622962 2579373 2254785 1982414 1794017 1626601 1502957 1333402 1100403 885109 714669 579468],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
  end;
  if t==1999 then
    plot2d(-1000*[48520 4290 3150 5330 9090 13130 16050 17560 17480 16380 16370 16430 17310]./[2462810 2596515 2612507 2315400 2014825 1777368 1555514 1403020 1235064 1004175 790910 624929 488707],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);  
    plot2d(1000*[41280 2840 1760 4100 11330 14910 12010 10440 9220 8700 9350 10600 12310]./[2453278 2598388 2614692 2311824 2009708 1815112 1637397 1513832 1366775 1135150 919534 730693 601520],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(-1000*[19787 1891 1649 4351 8632 13863 16252 16400 15159 14933 13830 14031 12660]./[2462810 2596515 2612507 2315400 2014825 1777368 1555514 1403020 1235064 1004175 790910 624929 488707],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);  
    plot2d(1000*[18084 1504 1305 3311 8227 12500 12172 10741 8866 8483 7736 8661 10027]./[2453278 2598388 2614692 2311824 2009708 1815112 1637397 1513832 1366775 1135150 919534 730693 601520],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
  end;
  if t==2000 then
    plot2d(-1000*[48260 4590 3120 5440 10250 16320 20670 21560 21130 18850 17920 17150 17800]./[2445814 2560634 2633243 2387921 2061216 1809314 1572366 1411341 1255354 1032921 813997 639575 501458],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);  
    plot2d(1000*[41500 3190 1770 4550 14680 20630 16620 13320 11150 9820 9980 10980 12630]./[2432565 2564711 2633417 2380628 2044320 1830479 1644042 1519941 1391941 1169661 952659 750344 620155],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(-1000*[20361 1996 1711 4301 8834 14995 18409 18474 17059 16036 15216 13871 14193]./[2445814 2560634 2633243 2387921 2061216 1809314 1572366 1411341 1255354 1032921 813997 639575 501458],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);  
    plot2d(1000*[18429 1592 1328 3454 9770 15549 15677 13481 10919 9510 9062 8841 11228]./[2432565 2564711 2633417 2380628 2044320 1830479 1644042 1519941 1391941 1169661 952659 750344 620155],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
  end;
  if t==2001 then
    plot2d(-1000*[47900 5100 3110 5380 11160 18890 25550 25230 24740 21350 19570 17940 18260]./[2435570 2518204 2634218 2474130 2121707 1843511 1591112 1416035 1267973 1058890 836525 654309 512382],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(1000*[41470 3750 1790 4630 17790 26300 21760 16070 13070 11010 10580 11380 12940]./[2419971 2523200 2635096 2457295 2088984 1841486 1646587 1520267 1407704 1203657 983773 774000 634731],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(-1000*[21324 2114 1743 4466 8915 16791 20835 20999 19242 17828 16833 14500 15059]./[2435570 2518204 2634218 2474130 2121707 1843511 1591112 1416035 1267973 1058890 836525 654309 512382],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
    plot2d(1000*[19349 1703 1457 3885 10821 19101 18570 15741 12777 10880 10091 9101 12046]./[2419971 2523200 2635096 2457295 2088984 1841486 1646587 1520267 1407704 1203657 983773 774000 634731],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
  end;
  if t==2002 then
    plot2d(-1000*[48500 5890 3170 5310 11900 20720 30290 28220 27970 23770 21240 18930 18820]./[2428333 2476101 2627409 2562592 2199201 1884553 1614076 1417338 1273583 1081520 858069 668846 521722],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(1000*[41970 4520 1820 4620 21240 32580 28720 19630 15570 12630 11410 11910 13290]./[2411137 2480728 2628053 2531106 2145881 1851250 1644165 1513904 1414006 1235712 1012918 800611 646169],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(-1000*[24150 2392 1860 4723 9549 18563 23769 23961 21435 19191 18523 15343 16118]./[2428333 2476101 2627409 2562592 2199201 1884553 1614076 1417338 1273583 1081520 858069 668846 521722],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
    plot2d(1000*[21830 1954 1480 4263 12377 23074 23296 19277 15378 12593 11181 9966 12658]./[2411137 2480728 2628053 2531106 2145881 1851250 1644165 1513904 1414006 1235712 1012918 800611 646169],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
  end;
  if t==2003 then
    plot2d(-1000*[49030 6880 3250 5240 12640 22430 35650 31620 31460 26520 23060 20080 19390]./[2419542 2441156 2610493 2645647 2295877 1936057 1642504 1417747 1273531 1099343 878056 683008 530115],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(1000*[42400 5500 1890 4360 23670 37650 35650 23020 17880 14210 12220 12450 13630]./[2401414 2443186 2612049 2596656 2217043 1863988 1636358 1500611 1411868 1263177 1040267 828667 656502],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(-1000*[26998 2765 1991 4818 10283 19896 27308 26273 24584 21921 20471 17107 17310]./[2419542 2441156 2610493 2645647 2295877 1936057 1642504 1417747 1273531 1099343 878056 683008 530115],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
    plot2d(1000*[24222 2189 1637 4523 14050 25950 27856 22446 18239 14363 12800 10930 13242]./[2401414 2443186 2612049 2596656 2217043 1863988 1636358 1500611 1411868 1263177 1040267 828667 656502],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
   end;
  if t==2004 then
    plot2d(-1000*[49600 7930 3420 5240 13050 23360 39330 33960 33390 28250 24280 21130 19880]./[2410082 2413955 2581959 2715937 2407919 1999593 1674022 1418946 1269245 1111103 896179 697017 538307],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(1000*[42870 6560 2060 4310 25150 41530 41900 26650 19910 15680 13080 12980 13970]./[2391299 2411635 2586487 2649984 2299904 1884111 1622895 1480943 1402579 1283256 1066489 856407 668299],[2.5:5:62.5],style=-9,rect=rrect,nax=nnax);
    plot2d(-1000*[29762 3162 2118 4652 10292 19654 28197 27987 26236 22901 20952 17908 16851]./[2410082 2413955 2581959 2715937 2407919 1999593 1674022 1418946 1269245 1111103 896179 697017 538307],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
    plot2d(1000*[26578 2784 1757 4573 14887 27210 30272 24847 20320 16106 13980 11946 13301]./[2391299 2411635 2586487 2649984 2299904 1884111 1622895 1480943 1402579 1283256 1066489 856407 668299],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
  end;
  if t==2005 then
    plot2d(-1000*[32003 3324 2133 4718 10373 19120 28539 29140 27198 24225 21325 19526 16705]./[2400652 2393233 2544612 2764641 2535108 2077890 1709531 1423408 1262752 1117155 912480 710936 546839],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
    plot2d(1000*[29043 2774 1839 4500 14671 26907 30897 25962 21256 17198 14831 13205 13162]./[2381255 2387101 2551350 2683348 2396257 1917306 1606211 1456503 1387651 1294646 1092013 882700 683323],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
  end;
  if t==2006 then
    plot2d(-1000*[33477 2999 2362 4806 10760 18827 28576 29177 27808 24882 22568 20458 16911]./[2391791 2379824 2500497 2787694 2670247 2169788 1750798 1434080 1255772 1118538 927665 725210 556028],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
    plot2d(1000*[29288 2528 1897 4552 14626 25853 30679 25831 21647 17787 15466 14059 13231]./[2371608 2370950 2507847 2696534 2498156 1966303 1591317 1430113 1367942 1297307 1117111 907341 702169],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
  end;
  if t==2007 then
    plot2d(-1000*[24638+7768 2854 2233 4859 10875 18405 28245 29258 26973 24761 22790 21316 17410]./[2382603 2369474 2455571 2795844 2798221 2275008 1800575 1452513 1248776 1116247 941623 739738 565776],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
    plot2d(1000*[21498+6968 2489 1892 4164 13664 24430 28959 24756 21108 17805 15567 14548 13422]./[2361866 2359037 2462683 2697479 2592647 2032353 1583299 1403260 1344096 1292124 1140919 930731 724095],[2.5:5:62.5],style=-5,rect=rrect,nax=nnax);
end;
endfunction
//************** HIV data *********************
function zz=showHIV()
  zz=0; clf(); xset("thickness",4); xset("font size",6); 
  for i=1:14,  
    x=[1:5]+5*(i-1);
    PrevHIVm(i)=sum(TIm(x)+TRm(x))/sum(Pm(x));
    PrevHIVf(i)=sum(TIf(x)+TRf(x))/sum(Pf(x));
  end;
  plot2d2(100*[-[PrevHIVm(1); PrevHIVm] [PrevHIVf(1); PrevHIVf]],[0:5:70]'*ones(1,2),style=[1 1],rect=[-50,0,50,60],nax=[0,5,0,7]);
 // plot2d2(100*[-(TIm+TRm)./Pm (TIf+TRf)./Pf],[0:99]'*ones(1,2),style=[1 1],rect=[-50,0,50,60],nax=[0,5,0,7]);
  
  plot2d([0 0],[0 xmax],style=1,rect=[-50,0,50,60],nax=[0,5,0,7]);
  xstring(30,55,string(t)); xstring(-40,10,"MEN"); xstring(20,10,"WOMEN"); xstring(46,2,"%"); xstring(-49,57,"AGE");
  if t==2002 then
    plot2d(-[5 4 8 22 24 18 12 12 5],[7.5 17.5:5:52.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]);
    plot2d(-[3 7; 2 6; 5 13; 15 32; 17 33; 11 28; 7 20; 7 20; 2 10]',ones(2,1)*[7.5 17.5:5:52.5],style=ones(1,9),rect=[-50,0,50,60],nax=[0,5,0,7]);
    plot2d([6 7.3 17.1 32.0 24.1 13.8 19.0 11.2 8],[7.5 17.5:5:52.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]);
    plot2d([3 9; 4.7 11.3; 12.9 22.3; 24.8 40.1; 17.3 32.5; 8.7 21.1; 12.8 27.2; 6.5 18.7; 4 15]',ones(2,1)*[7.5 17.5:5:52.5],style=ones(1,9),rect=[-50,0,50,60],nax=[0,5,0,7]);
  end;
  if t==2005 then 
    plot2d(-[4.9 4.2 1.6 3.2 6.0 12.1 23.3 23.3 17.5 10.3 14.2 6.4],[2.5:5:57.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]); 
    plot2d(-[1.8 12.8; 2.2 8.0; 0.8 3.4; 1.4 7.1; 3.8 9.4; 8.0 17.9; 17.2 30.7; 17.8 29.8; 12.0 24.7; 6.5 16.0; 8.5 22.7; 1.9 19.8]',ones(2,1)*[2.5:5:57.5],style=ones(1,12),rect=[-50,0,50,60],nax=[0,5,0,7]);
    plot2d([5.3 4.8 1.8 9.4 23.9 33.3 26.0 19.3 12.4 8.7 7.5 3.0],[2.5:5:57.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]); 
    plot2d([3.1 9.0; 2.8 8.0; 0.9 3.4; 7.1 12.4; 19.8 28.4; 27.7 39.4; 21.5 30.9; 14.9 24.6; 9.4 16.2; 6.0 12.6; 4.9 11.2; 1.6 5.6]',ones(2,1)*[2.5:5:57.5],style=ones(1,12),rect=[-50,0,50,60],nax=[0,5,0,7]);
  end;
  if t==2008 then
    plot2d(-[3.0 2.5 5.1 15.7 25.8 18.5 19.2 8.4 10.4 6.2],[7.5 17.5:5:57.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]);
    plot2d(-[2 5; 1.1 7.5; 3.5 8; 11 22; 19 34; 13 26; 14 27; 5 13; 6 16; 3.5 11]',ones(2,1)*[7.5 17.5:5:57.5],style=ones(1,10),rect=[-50,0,50,60],nax=[0,5,0,7]);
    plot2d([2.0 6.7 21.1 32.7 29.1 24.8 16.3 14.1 10.2 7.7],[7.5 17.5:5:57.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]);
    plot2d([1 3; 5 9; 17.5 25.5; 27 38; 24 36; 20 30; 13 21; 10 20; 7 15.5; 5 12.5]',ones(2,1)*[7.5 17.5:5:57.5],style=ones(1,10),rect=[-50,0,50,60],nax=[0,5,0,7]);
  end;
  if t==1991 then plot2d([2 1.8 1.2 0.6 0.3 0.6 0],[17.5:5:47.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]); end;
  if t==1992 then plot2d([3 3.7 2.5 2.2 1.8 2.2 0],[17.5:5:47.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]); end;
  if t==1993 then plot2d([5 6.2 4.9 3.1 1.5 2.2 0],[17.5:5:47.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]); end;
  if t==1994 then plot2d([7 9.2 8.9 6.5 3.7 5.2 0],[17.5:5:47.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]); end;
  if t==1995 then plot2d([10 12.9 10.8 8 6.8 4.3 7.5],[17.5:5:47.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]); end;
  if t==1996 then plot2d([13 17.8 15.1 12.3 9.2 9.8 6],[17.5:5:47.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]); end;
  if t==1997 then plot2d([13 20 17.8 14.2 9.2 7.4 9],[17.5:5:47.5],style=-9,rect=[-50,0,50,60],nax=[0,5,0,7]); end;
  if t==1998 then
    plot2d([21 26.1 26.9 19.1 13.4 10.5 10.2],[17.5:5:47.5],style=-9);
    plot2d([18.4 23.8; 24.1 28.1; 24.7 29.0; 17.1 21.1; 11.2 15.6; 6.8 14.1; 0.4 20.0]',ones(2,1)*[17.5:5:47.5],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==1999 then
    plot2d([16.5 25.6 26.4 21.7 16.2 12.0 7.5],[17.5:5:47.5],style=-9);
    plot2d([14.9 18.1; 24.0 27.3; 24.6 28.3; 19.1 23.8; 14.1 18.3; 8.5 15.6; 0 15.9]',ones(2,1)*[17.5:5:47.5],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==2000 then
    plot2d([16.1 29.1 30.6 23.3 15.8 10.2 13.1],[17.5:5:47.5],style=-9);
    plot2d([14.5 17.7; 27.4 30.8; 28.8 32.4; 21.5 25.1; 13.9 17.7; 6.9 13.3; 2.09 24.0]',ones(2,1)*[17.5:5:47.5],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==2001 then
    plot2d([15.4 28.4 31.4 25.6 19.3 9.8],[17.5:5:37.5 45],style=-9);
    plot2d([13.8 16.9; 26.5 30.2; 29.5 33.3; 23.5 27.7; 17.0 21.5; 7.0 12.6]',ones(2,1)*[17.5:5:37.5 45],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==2002 then
    plot2d([14.8 29.1 34.5 29.5 19.8 17.2],[17.5:5:37.5 45],style=-9);
    plot2d([13.4 16.1; 27.5 30.6; 32.6 36.4; 27.4 31.6; 17.5 22.0; 13.5 20.9]',ones(2,1)*[17.5:5:37.5 45],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==2003 then
    plot2d([15.8 30.3 35.4 30.9 23.4 15.8],[17.5:5:37.5 45],style=-9);
    plot2d([14.3 17.2; 28.8 31.8; 33.6 37.2; 28.9 32.9; 20.9 25.9; 12.3 19.3]',ones(2,1)*[17.5:5:37.5 45],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==2004 then
    plot2d([16.1 30.8 38.5 34.4 24.5 17.5],[17.5:5:37.5 45],style=-9);
    plot2d([14.7 17.5; 29.3 32.3; 36.8 40.3; 32.2 36.6; 21.9 27.2; 14.0 21.0]',ones(2,1)*[17.5:5:37.5 45],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==2005 then
    plot2d([15.9 30.6 39.5 36.4 28.0 19.8],[17.5:5:37.5 45],style=-9);
    plot2d([14.6 17.2; 29.0 32.2; 37.7 41.3; 34.3 38.5; 25.2 30.8; 16.1 23.6]',ones(2,1)*[17.5:5:37.5 45],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==2006 then
    plot2d([13.7 28.0 38.7 37.0 29.3 21.3 15.5],[17.5:5:47.5],style=-9);
    plot2d([12.8 14.6; 26.9 29.1; 37.3 40.2; 35.5 38.5; 27.7 31.5; 18.5 24.1; 8.5 25.0]',ones(2,1)*[17.5:5:47.5],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==2007 then
    plot2d([13.1 28.0 37.5 39.6 33.0 22.2 20.6],[17.5:5:47.5],style=-9);
    plot2d([12.2 14.0; 26.9 29.1; 36.2 38.8; 38.0 41.2; 31.1 34.9; 19.1 25.7; 13.2 30.7]',ones(2,1)*[17.5:5:47.5],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
  end;
  if t==2008 then
    plot2d([14.1 26.9 37.9 40.4 32.4 23.3 17.6],[17.5:5:47.5],style=-9);
    plot2d([13.1 15.0; 25.9 27.9; 36.4 39.3; 38.7 42.0; 30.5 34.3; 20.3 26.6; 10.7 27.7]',ones(2,1)*[17.5:5:47.5],style=ones(1,7),rect=[-50,0,50,60],nax=[0,5,0,7]); 
end;
endfunction
//**************** show assumptions *******************
if showassumptions==1 then
  clf(); xset("thickness",3); xset("font size",4);
  chdir(loadDirectory);
  dataage=read('lifetime.txt',21,17);
  moyenne=zeros(xmax,1);
  for x=15:30,
    moyenne(x)=sum(dataage(2:21,1).*dataage(2:21,x-13))/sum(dataage(2:21,x-13));  
  end;
  plot2d([15:30],moyenne(15:30),style=-9,rect=[0,0,60,10],nax=[0,7,0,3]);
  plot2d(1:xmax,cumsum(rate/underreport),rect=[0,0,60,10],nax=[0,7,0,3]);
  xstring(55,0.1,"AGE");
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"cumul.eps"); end; xclick(); 
  clf(); xset("thickness",3); xset("font size",4);
  plot2d([1:xmax],rate/underreport,rect=[0,0,60,1],nax=[0,7,1,2]);
  xstring(1,0.95,"/YR"); xstring(55,0.1,"AGE");
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"rate.eps"); end; xclick(); 
  AIDSsurvival0f(1)=1; AIDSsurvival20f(1)=1; 
  for y=2:25,
    AIDSsurvival0f(y)=(1-df(y-1))*(1-AIDS(y-1,y-1))*AIDSsurvival0f(y-1);
    AIDSsurvival20f(y)=(1-df(20+y-1))*(1-AIDS(20+y-1,y-1))*AIDSsurvival20f(y-1);
  end;
  clf(); xset("thickness",3); xset("font size",4);
  plot2d([0:24]'*ones(1,2),100*[AIDSsurvival0f AIDSsurvival20f],style=[1 1],rect=[0,0,20,100],nax=[0,3,0,3]); 
  xstring(0.5,93,"%"); xstring(17,9,"YEARS"); 
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"AIDSsurvival.eps"); end; xclick();
  clf(); xset("thickness",3); xset("font size",4);
  plot2d([1:xmax],1000*dm,rect=[0,0,80,50],nax=[0,5,0,6]);
  plot2d([1:xmax],1000*df,rect=[0,0,80,50],nax=[0,5,0,6]);
  xstring(2,47,"/1000/YEAR"); xstring(73,1,"AGE"); 
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"mortality.eps"); end; xclick();
  for t=t0:tmax,
    z=condom(t);
    con20(t-t0+1)=z(20);
    con35(t-t0+1)=z(35);
    con60(t-t0+1)=z(60);
    mother(t-t0+1)=M2C(t);
  end;
  clf(); xset("thickness",3); xset("font size",4);
  plot2d([t0:tmax],100*[con20 con35 con60],style=[1 1 1]);
  condo1524_1998=((0.186*28+0.212*827+0.195*854)/(28+827+854)+(0.093*274+0.187*1354+0.144*1628)/(274+1354+1628))/2;
  plot2d([1998 2002 2005 2008],[100*condo1524_1998 46.1 55.7 73.1],style=-9);
  plot2d(ones(2,1)*[2002 2005 2008],[42 50; 51 60; 69 77]',style=ones(1,3));
  condo2549_1998=(0.088*660+0.146*937+0.076*1597+0.065*805+0.171*627+0.066*1433+0.069*910+0.090*481+0.026*1391+0.043*678+0.107*339+0.035*1018+0.026*511+0.14*185+0.03*696)/(660+937+1597+805+627+1433+910+481+1391+678+339+1018+511+185+696);
  plot2d([1998 2002 2005 2008],[100*condo2549_1998 19.7 29.1 58.1],style=-4,rect=[1985,0,tmax,100],nax=[0,8,0,3]);
  plot2d(ones(2,1)*[2002 2005 2008],[18 22; 28 32; 55 61]',style=ones(1,3));
  plot2d([2002 2005 2008],[5.6 5.3 25.9],style=-5,rect=[1985,0,tmax,100],nax=[0,8,0,3]);
  plot2d(ones(2,1)*[2002 2005 2008],[3 9; 3 9; 20 35]',style=ones(1,3),rect=[1985,0,tmax,100],nax=[0,8,0,3]);
  xstring(1986,95,"%"); 
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"condom.eps"); end; xclick();
  clf(); xset("thickness",3); xset("font size",4);
  plot2d([20 20 20 20],[100*condo1524_1998 46.1 55.7 73.1],style=-9,rect=[0,0,80,100],nax=[0,5,0,3]);
  plot2d([35 35 35 35],[100*condo2549_1998 19.7 29.1 58.1],style=-4,rect=[0,0,80,100],nax=[0,5,0,3]);
  plot2d([60 60 60],[5.6 5.3 25.9],style=-5,rect=[0,0,80,100],nax=[0,5,0,3]);
  c1998=condom(1998); c2002=condom(2002); c2005=condom(2005); c2008=condom(2008); 
  plot2d([debut:xmax]'*ones(1,4),100*[c1998(debut:xmax) c2002(debut:xmax) c2005(debut:xmax) c2008(debut:xmax)],style=ones(1,4),rect=[0,0,80,100],nax=[0,5,0,3]);
  plot2d([20 20; 20 20; 20 20; 35 35; 35 35; 35 35; 60 60; 60 60; 60 60]',[42 50; 51 60; 69 77; 18 22; 28 32; 55 61; 3 9; 3 9; 20 35]',style=ones(1,9),rect=[0,0,80,100],nax=[0,5,0,3]);
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"condom2.eps"); end; xclick();  
  clf(); xset("thickness",3); xset("font size",4);
  plot2d2([1:xmax],100*weight/sum(weight),rect=[0,0,60,6],nax=[0,4,0,0]);
  xstring(52,0.1,"AGE"); 
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"fertility.eps"); end; xclick();
  clf(); xset("thickness",3); xset("font size",4);
  plot2d([debut:xmax]'*ones(1,3),[ym(debut:xmax)' yave(debut:xmax)' yM(debut:xmax)'],style=[1 1 1],rect=[0,0,80,80],nax=[0,5,0,5]);
  chdir(loadDirectory);
  data=read('couples.txt',-1,2);
  plot2d(data(:,1),data(:,2),style=-1,rect=[0,0,80,80],nax=[0,5,0,5]);
  xstring(55,1,"WOMEN''S AGE"); xstring(1,75,"MEN''S AGE"); 
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"gug.eps"); end; xclick();
  clf(); xset("thickness",3); xset("font size",4);
  plot2d([debut:xmax]'*ones(1,3),[ym(debut:xmax)' yave(debut:xmax)' yM(debut:xmax)'],style=[1 1 1],rect=[0,0,80,80],nax=[0,5,0,5]);
  chdir(loadDirectory);
  data2=read('couples_2.txt',-1,2);
  plot2d(data2(:,1),data2(:,2),style=-1,rect=[0,0,80,80],nax=[0,5,0,5]);
  xstring(55,1,"WOMEN''S AGE"); xstring(1,75,"MEN''S AGE"); 
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"carltonville.eps"); end; xclick();
  clf();  xset("thickness",3); xset("font size",4);
  plot2d([-Sm0(1:5:96) Sf0(1:5:96)],[2:5:97]'*ones(1,2),style=[-9 -9],rect=[-25,0,25,80],nax=[0,3,0,5]);  
  plot2d2([-Sm Sf],[0:99]'*ones(1,2),style=[1 1],rect=[-25,0,25,80],nax=[0,3,0,5]);   
  plot2d([0 0],[0 80],rect=[-25,0,25,80],nax=[0,3,0,5]); xstring(-24,75,"AGE"); xstring(-20,50,"MEN"); xstring(10,50,"WOMEN");
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"agepyramid.eps"); end; xclick();
  clf();  xset("thickness",3); xset("font size",4);
  plot2d([2003 2008 tmax],[0 75 75],style=1,rect=[t0,0,tmax,100],nax=[0,8,0,3]);
  plot2d([2004:2008],[15 34 50 57 75],style=-9); xstring(t0+1,95,"%");
  plot2d(ones(2,1)*[2004:2008],[13 17; 29 40; 43 60; 49 69; 55 95]',style=ones(1,5),rect=[1985,0,tmax,100],nax=[0,8,0,3]); 
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"M2C.eps"); end; xclick();
  clf(); xset("thickness",3); xset("font size",4);
  xset("line style",2);
  plot2d([2003 2007 tmax],[0 100 100],style=1,rect=[2000,0,tmax,100],nax=[0,5,0,3]);
  xset("line style",1);
  xstring(2000.5,95,"%");
  plot2d([2010 2012 tmax],[0 100 100],style=1,rect=[2000,0,tmax,100],nax=[0,5,0,3]); 
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"ARTaccess.eps"); end; xclick();
  clf(); xset("thickness",3); xset("font size",4); 
  plot2d([t0:2020],liste([t0:2020]-t0+1)/liste(1),rect=[1985,0,2020,1],nax=[0,8,0,2]);
  if savefig==1 then chdir(saveDirectory); xs2eps(0,"births.eps"); end; xclick();
end;
//**************** program **************************
If(16,1)=1;
Sf(16)=Sf(16)-1;
for t=t0:tmax,
  TIf=sum(If,'c'); TRf=sum(Rf,'c'); Pf=Sf+TIf+TRf;
  TIm=sum(Im,'c'); TRm=sum(Rm,'c'); Pm=Sm+TIm+TRm;
  //*** output ***
  TotalPopulation(t-t0+1)=sum(Pm+Pf);
  TotalART(t-t0+1)=sum(TRm+TRf)/sum(Pm+Pf);
  CrudeDeathRate(t-t0+1)=(sum(dm.*(Sm+TIm+TRm))+sum((Im+If).*AIDS)+sum((Rm+Rf).*ARTdeath)+sum(df.*(Sf+TIf+TRf)))/sum(Pm+Pf);
  incidence(t-t0+1)=sum(Im(14:99,1)+If(14:99,1))/sum(Pm(14:99)+Pf(14:99));
  incidence1524(t-t0+1)=sum(Im(16:25,1)+If(16:25,1))/sum(Pm(16:25)+Pf(16:25));
  incidence2534(t-t0+1)=sum(Im(26:35,1)+If(26:35,1))/sum(Pm(26:35)+Pf(26:35));
  incidence3544(t-t0+1)=sum(Im(36:45,1)+If(36:45,1))/sum(Pm(36:45)+Pf(36:45));
  incidence4554(t-t0+1)=sum(Im(46:55,1)+If(46:55,1))/sum(Pm(46:55)+Pf(46:55));
  incidence5564(t-t0+1)=sum(Im(56:65,1)+If(56:65,1))/sum(Pm(56:65)+Pf(56:65));
  HIV1549(t-t0+1)=sum(TIm(16:50)+TRm(16:50)+TIf(16:50)+TRf(16:50))/sum(Pf(16:50)+Pm(16:50));
  HIVnoART1549(t-t0+1)=sum(TIm(16:50)+TIf(16:50))/sum(Pf(16:50)+Pm(16:50));
  if show=='HIV' then zz=showHIV(); end;
  if show=='PYRAMID' then zz=showpyramid(); end;
  if show=='DEATH' then zz=showdeath(); end;
  if savefig==1 then chdir(saveDirectory); xs2eps(0,show+string(t)+".eps"); end;
  //*** end output ***
  infected_mother=sum(weight.*TIf)/(sum(weight.*Pf));
  safe_child=1-M2C(t)*infected_mother; 
  SSm(1)=B(t)*safe_child;
  SSf(1)=B(t)*safe_child;
  IIm(1,1)=B(t)*(1-safe_child);
  IIf(1,1)=B(t)*(1-safe_child);
  ratecondom=(1-condom(t)).*rate;
  HAART=ART(t);
  for x=1:xmax-1,
    incidencef(x)=1-exp(-pf*ratecondom(x)*choix(x,:)*((TIm+epsi*TRm)./Pm));
    incidencem(x)=1-exp(-pm*choix(:,x)'*(ratecondom.*(TIf+epsi*TRf))/Pm(x));
    SSm(x+1)=(1-dm(x))*Sm(x)*(1-incidencem(x));
    SSf(x+1)=(1-df(x))*Sf(x)*(1-incidencef(x));
    IIm(x+1,1)=(1-dm(x))*Sm(x)*incidencem(x);
    IIf(x+1,1)=(1-df(x))*Sf(x)*incidencef(x);
    IIm(x+1,2:x+1)=(1-dm(x))*Im(x,1:x).*(1-AIDS(x,1:x)).*(1-HAART(x,1:x))+drop*Rm(x,1:x)*(1-dm(x)).*(1-ARTdeath(x,1:x));
    IIf(x+1,2:x+1)=(1-df(x))*If(x,1:x).*(1-AIDS(x,1:x)).*(1-HAART(x,1:x))+drop*Rf(x,1:x)*(1-df(x)).*(1-ARTdeath(x,1:x));
    RRm(x+1,1:x)=(1-drop)*Rm(x,1:x)*(1-dm(x)).*(1-ARTdeath(x,1:x))+Im(x,1:x)*(1-dm(x)).*(1-AIDS(x,1:x)).*HAART(x,1:x);
    RRf(x+1,1:x)=(1-drop)*Rf(x,1:x)*(1-df(x)).*(1-ARTdeath(x,1:x))+If(x,1:x)*(1-df(x)).*(1-AIDS(x,1:x)).*HAART(x,1:x);
  end;
  if t==1990 then inc_f_1990=incidencef; inc_m_1990=incidencem; end;
  if t==2000 then inc_f_2000=incidencef; inc_m_2000=incidencem; end;
  if t==2010 then inc_f_2010=incidencef; inc_m_2010=incidencem; end;
  Sm=SSm; Sf=SSf; Im=IIm; If=IIf; Rm=RRm; Rf=RRf;
end;
//***************** show aggregated results ************************
clf(); xset("thickness",3); xset("font size",4);
plot2d([t0:tmax],TotalPopulation(1:tmax-t0+1)/TotalPopulation(1),rect=[1985,0,tmax,2],nax=[0,8,0,3]);
plot2d([1986:2008],[35101 35948 36795 37639 38476 39300 40091 40941 41662 42228 42794  43354 43962 44526 45064 45576 46077 46567 47033 47483 47926 48367 48783]/35101,style=-9,rect=[1985,0,tmax,2],nax=[0,8,0,3]);
//if savefig==1 then chdir(saveDirectory2); xs2eps(0,"TotalPopulation.eps"); end; 
xclick();
clf(); xset("thickness",3); xset("font size",4);
plot2d([t0:tmax],100*TotalART(1:tmax-t0+1),rect=[1990,0,tmax,10],nax=[1,4,0,3]);
plot2d([2004+10/12 2005+8/12 2006+6/12 2007+11/12 2008+11/12],100*[19500 78000 140000 458951 700500]./[47.032 47.483 47.925 48.367 48.783]/(1e6),style=-9,rect=[1990,0,tmax,10],nax=[1,4,0,3]);
xstring(1990.5,9.5,"%");
//if savefig==1 then xs2eps(0,"ARTdata.eps"); end; 
xclick();
clf(); xset("thickness",3); xset("font size",4);
plot2d([t0:tmax],1000*CrudeDeathRate(1:tmax-t0+1),rect=[1990,0,tmax,20]);
plot2d([1997:2004],1000*([217630 235620 253470 276620 289990 321020 345110 365800]+[176980 192170 209120 232550 255790 284890 311370 336620])./[43.354 43.962 44.526 45.064 45.576 46.077 46.567 47.033]/(1e6),style=-9,rect=[1990,0,tmax,20]); 
plot2d([1997:2007],1000*[317131 365852 381820 415983 454847 502031 556769 576700 598054 612462 601133]./[43.354 43.962 44.526 45.064 45.576 46.077 46.567 47.033 47.483 47.926 48.367]/(1e6),style=-5,rect=[1990,0,tmax,20]); 
plot2d([1985:2007],[8.9 8.6 8.26 7.96 7.70 7.5 7.34 7.26 7.26 7.36 7.61 8.05 8.68 9.51 10.51 11.64 12.82 14.0 15.09 15.96 16.42 16.57 16.66],style=-4,rect=[1990,0,tmax,20],nax=[1,4,0,5]);
xstring(1990.5,19,"/1000");
//if savefig==1 then xs2eps(0,"CrudeDeathData.eps"); end; 
xclick();
clf(); xset("thickness",3); xset("font size",4);
plot2d([t0:tmax],100*HIV1549(1:tmax-t0+1),rect=[1985,0,tmax,40],nax=[0,8,0,5]);
xset("line style",1);
plot2d([t0:tmax],100*HIVnoART1549(1:tmax-t0+1),rect=[1985,0,tmax,40],nax=[0,8,0,5]);
xset("line style",1);
plot2d([1990:2008],...//(16.9/29.3)*
[0.7 1.7 2.2 4.0 7.6 10.4 14.2 17.0 22.8 22.4 24.5 24.8 26.5 27.9 29.5 30.2 29.1 29.4 29.3],style=-9,rect=[1985,0,tmax,40]);
//plot2d([2008]*ones(1,2),[28.5 30.1],rect=[1985,0,tmax,40],nax=[0,8,0,5]); 
plot2d([2002 2005 2008],[15.6 16.2 16.9],style=-4,rect=[1985,0,tmax,40],nax=[0,8,0,5]); 
plot2d(ones(2,1)*[2002 2005 2008],[13.9 17.6; 14.9 17.7; 15.5 18.4]',style=ones(1,3),rect=[1985,0,tmax,40],nax=[0,8,0,5]); 
xstring(1986,38,"%");
//if savefig==1 then xs2eps(0,"ANC.eps"); end; 
xclick();
clf(); xset("thickness",3); xset("font size",4);
plot2d([t0:tmax],100*incidence,rect=[1990,0,2020,2],nax=[0,4,0,3]);
xstring(1990.5,1.85,"%/YR");
plot2d(2008,100*380000/35000000,style=-9);
plot2d(2008,100*430000/35000000,style=-4,rect=[1990,0,2020,2],nax=[0,4,0,3]);
//plot2d([2008 2008],100*[350000 390000]/35000000,rect=[1990,0,2020,2],nax=[0,4,0,3]);
//if savefig==1 then xs2eps(0,"incidence.eps"); end; 
xclick();
//clf(); xset("thickness",3); xset("font size",4);
//plot2d([t0:tmax]',100*[incidence1524 incidence2534 incidence3544 incidence4554 incidence5564],style=[1 2 5 3 4],rect=[1990,0,2020,4],nax=[0,4,0,3]);
//xstring(1990.5,3.8,"%/YR"); xclick();
clf(); xset("thickness",3); xset("font size",4);
plot2d([1:xmax-1]',100*[inc_f_1990 inc_f_2000 inc_f_2010],style=[1 1 1],rect=[15,0,65,6],nax=[0,6,0,4]);
xstring(20,0.7,"1990"); xstring(20,2.7,"2010"); xstring(20,5.2,"2000"); xstring(15.5,5.7,"%/YR"); xstring(40,5,"WOMEN");
xclick();
clf(); xset("thickness",3); xset("font size",4);
plot2d([1:xmax-1]',100*[inc_m_1990 inc_m_2000 inc_m_2010],style=[1 1 1],rect=[15,0,65,6],nax=[0,6,0,4]);
xstring(30,0.35,"1990"); xstring(30,3.3,"2000"); xstring(30,1.9,"2010"); xstring(15.5,5.7,"%/YR"); xstring(40,5,"MEN");
xclick();
clf(); xset("thickness",3); xset("font size",4);
plot2d([17.5:5:67.5],[13.5 20.5 24 17 12 7 3.5 2.35 1.2 0 0],rect=[15,0,65,30],nax=[0,6,0,0]);
plot2d([17.5:5:67.5],[30 25 19 11 6 4 2 1 0.5 0 0],rect=[15,0,65,30],nax=[0,6,0,0]);
xstring(60,0.4,"AGE"); xstring(25,8,"WOMEN"); xstring(35,15,"MEN");