--- imach/src/imach.c	2015/07/16 16:49:02	1.192
+++ imach/src/imach.c	2015/09/15 17:34:58	1.201
@@ -1,6 +1,43 @@
-/* $Id: imach.c,v 1.192 2015/07/16 16:49:02 brouard Exp $
+/* $Id: imach.c,v 1.201 2015/09/15 17:34:58 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.201  2015/09/15 17:34:58  brouard
+  Summary: 0.98r0
+
+  - Some new graphs like suvival functions
+  - Some bugs fixed like model=1+age+V2.
+
+  Revision 1.200  2015/09/09 16:53:55  brouard
+  Summary: Big bug thanks to Flavia
+
+  Even model=1+age+V2. did not work anymore
+
+  Revision 1.199  2015/09/07 14:09:23  brouard
+  Summary: 0.98q6 changing default small png format for graph to vectorized svg.
+
+  Revision 1.198  2015/09/03 07:14:39  brouard
+  Summary: 0.98q5 Flavia
+
+  Revision 1.197  2015/09/01 18:24:39  brouard
+  *** empty log message ***
+
+  Revision 1.196  2015/08/18 23:17:52  brouard
+  Summary: 0.98q5
+
+  Revision 1.195  2015/08/18 16:28:39  brouard
+  Summary: Adding a hack for testing purpose
+
+  After reading the title, ftol and model lines, if the comment line has
+  a q, starting with #q, the answer at the end of the run is quit. It
+  permits to run test files in batch with ctest. The former workaround was
+  $ echo q | imach foo.imach
+
+  Revision 1.194  2015/08/18 13:32:00  brouard
+  Summary:  Adding error when the covariance matrix doesn't contain the exact number of lines required by the model line.
+
+  Revision 1.193  2015/08/04 07:17:42  brouard
+  Summary: 0.98q4
+
   Revision 1.192  2015/07/16 16:49:02  brouard
   Summary: Fixing some outputs
 
@@ -677,11 +714,12 @@ typedef struct {
 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
 #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
 #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
-#define codtabm(h,k)  1 & (h-1) >> (k-1) ;
+#define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
 #define MAXN 20000
 #define YEARM 12. /**< Number of months per year */
 #define AGESUP 130
 #define AGEBASE 40
+#define AGEOVERFLOW 1.e20
 #define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */
 #ifdef _WIN32
 #define DIRSEPARATOR '\\'
@@ -693,11 +731,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.192 2015/07/16 16:49:02 brouard Exp $ */
+/* $Id: imach.c,v 1.201 2015/09/15 17:34:58 brouard Exp $ */
 /* $State: Exp $ */
-
-char version[]="Imach version 0.98q3, July 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
-char fullversion[]="$Revision: 1.192 $ $Date: 2015/07/16 16:49:02 $"; 
+#include "version.h"
+char version[]=__IMACH_VERSION__;
+char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
+char fullversion[]="$Revision: 1.201 $ $Date: 2015/09/15 17:34:58 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -763,7 +802,7 @@ char command[FILENAMELENGTH];
 int  outcmd=0;
 
 char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
-
+char fileresu[FILENAMELENGTH]; /* Without r in front */
 char filelog[FILENAMELENGTH]; /* Log file */
 char filerest[FILENAMELENGTH];
 char fileregp[FILENAMELENGTH];
@@ -831,7 +870,7 @@ int estepm;
 
 int m,nb;
 long *num;
-int firstpass=0, lastpass=4,*cod, *Tage,*cens;
+int firstpass=0, lastpass=4,*cod, *cens;
 int *ncodemax;  /* ncodemax[j]= Number of modalities of the j th
 		   covariate for which somebody answered excluding 
 		   undefined. Usually 2: 0 and 1. */
@@ -851,8 +890,9 @@ double  **covar; /**< covar[j,i], value
 		  * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
 double  idx; 
 int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
+int *Tage;
 int *Ndum; /** Freq of modality (tricode */
-int **codtab; /**< codtab=imatrix(1,100,1,10); */
+/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
 int **Tvard, *Tprod, cptcovprod, *Tvaraff;
 double *lsurv, *lpop, *tpop;
 
@@ -1888,13 +1928,16 @@ double **prevalim(double **prlim, int nl
     if(nagesqr==1)
       cov[3]= agefin*agefin;;
     for (k=1; k<=cptcovn;k++) {
-      cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
-      /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/
+      /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+      cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+      /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
     }
     /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-    for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2];
+    /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
+    for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
     for (k=1; k<=cptcovprod;k++) /* Useless */
-      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+      /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
     
     /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
     /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
@@ -2067,12 +2110,15 @@ double ***hpxij(double ***po, int nhstep
       if(nagesqr==1)
 	cov[3]= agexact*agexact;
       for (k=1; k<=cptcovn;k++) 
-	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+	/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
       for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
 	/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
+	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+	/* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
       for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
-	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+	/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
 
 
       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
@@ -2530,7 +2576,7 @@ void likelione(FILE *ficres,double p[],
   int k;
 
   if(*globpri !=0){ /* Just counts and sums, no printings */
-    strcpy(fileresilk,"ilk"); 
+    strcpy(fileresilk,"ILK_"); 
     strcat(fileresilk,fileres);
     if((ficresilk=fopen(fileresilk,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresilk);
@@ -2547,7 +2593,7 @@ void likelione(FILE *ficres,double p[],
   *fretone=(*funcone)(p);
   if(*globpri !=0){
     fclose(ficresilk);
-    fprintf(fichtm,"\n<br>File of contributions to the likelihood: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
+    fprintf(fichtm,"\n<br>File of contributions to the likelihood (if mle=1): <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
     fflush(fichtm); 
   } 
   return;
@@ -2580,7 +2626,7 @@ void mlikeli(FILE *ficres,double p[], in
     for (j=1;j<=npar;j++)
       xi[i][j]=(i==j ? 1.0 : 0.0);
   printf("Powell\n");  fprintf(ficlog,"Powell\n");
-  strcpy(filerespow,"pow"); 
+  strcpy(filerespow,"POW_"); 
   strcat(filerespow,fileres);
   if((ficrespow=fopen(filerespow,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", filerespow);
@@ -2891,7 +2937,7 @@ void lubksb(double **a, int n, int *indx
 
 void pstamp(FILE *fichier)
 {
-  fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart);
+  fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
 }
 
 /************ Frequencies ********************/
@@ -2907,8 +2953,8 @@ void  freqsummary(char fileres[], int ia
   
   pp=vector(1,nlstate);
   prop=matrix(1,nlstate,iagemin,iagemax+3);
-  strcpy(fileresp,"p");
-  strcat(fileresp,fileres);
+  strcpy(fileresp,"P_");
+  strcat(fileresp,fileresu);
   if((ficresp=fopen(fileresp,"w"))==NULL) {
     printf("Problem with prevalence resultfile: %s\n", fileresp);
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
@@ -2943,13 +2989,13 @@ void  freqsummary(char fileres[], int ia
 	bool=1;
 	if  (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
 	  for (z1=1; z1<=cptcoveff; z1++)       
-            if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]){
+            if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
                 /* Tests if the value of each of the covariates of i is equal to filter j1 */
               bool=0;
-              /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n", 
-                bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1],
-                j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);*/
-              /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/
+              /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
+                bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
+                j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
+              /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
             } 
 	}
  
@@ -2978,10 +3024,10 @@ void  freqsummary(char fileres[], int ia
       pstamp(ficresp);
       if  (cptcovn>0) {
 	fprintf(ficresp, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresp, "**********\n#");
 	fprintf(ficlog, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficlog, "**********\n#");
       }
       for(i=1; i<=nlstate;i++) 
@@ -3109,7 +3155,7 @@ void prevalence(double ***probs, double
 	bool=1;
 	if  (cptcovn>0) {
 	  for (z1=1; z1<=cptcoveff; z1++) 
-	    if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) 
+	    if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) 
 	      bool=0;
 	} 
 	if (bool==1) { 
@@ -3375,14 +3421,15 @@ void tricode(int *Tvar, int **nbcode, in
        nbcode[Tvar[j]][1]=0;
        nbcode[Tvar[j]][2]=1;
        nbcode[Tvar[j]][3]=2;
+       To be continued (not working yet).
     */
-    ij=0; /* ij is similar to i but can jumps over null modalities */
-    for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 to 1*/
-    	if (Ndum[i] == 0) { /* If at least one individual responded to this modality k */
+    ij=0; /* ij is similar to i but can jump over null modalities */
+    for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
+    	if (Ndum[i] == 0) { /* If nobody responded to this modality k */
 	  break;
 	}
 	ij++;
-	nbcode[Tvar[j]][ij]=i;  /* stores the original modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
+	nbcode[Tvar[j]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
 	cptcode = ij; /* New max modality for covar j */
     } /* end of loop on modality i=-1 to 1 or more */
       
@@ -3775,11 +3822,11 @@ void varevsij(char optionfilefiname[], d
 
   if(popbased==1){
     if(mobilav!=0)
-      strcpy(digitp,"-populbased-mobilav-");
-    else strcpy(digitp,"-populbased-nomobil-");
+      strcpy(digitp,"-POPULBASED-MOBILAV_");
+    else strcpy(digitp,"-POPULBASED-NOMOBIL_");
   }
   else 
-    strcpy(digitp,"-stablbased-");
+    strcpy(digitp,"-STABLBASED_");
 
   if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
@@ -3789,7 +3836,7 @@ void varevsij(char optionfilefiname[], d
     }
   }
 
-  strcpy(fileresprobmorprev,"prmorprev"); 
+  strcpy(fileresprobmorprev,"PRMORPREV-"); 
   sprintf(digit,"%-d",ij);
   /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
   strcat(fileresprobmorprev,digit); /* Tvar to be done */
@@ -3812,7 +3859,8 @@ void varevsij(char optionfilefiname[], d
   }  
   fprintf(ficresprobmorprev,"\n");
   fprintf(ficgp,"\n# Routine varevsij");
-  /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
+  fprintf(ficgp,"\nunset title \n");
+/* fprintf(fichtm, "#Local time at start: %s", strstart);*/
   fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
   fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
 /*   } */
@@ -4024,9 +4072,11 @@ void varevsij(char optionfilefiname[], d
   free_vector(gmp,nlstate+1,nlstate+ndeath);
   free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
   free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
-  fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240");
+  /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */
+  fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480");
   /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
   fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
+  fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
 /*   fprintf(ficgp,"\n plot \"%s\"  u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
@@ -4034,11 +4084,11 @@ void varevsij(char optionfilefiname[], d
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
   fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
-  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
-  /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);
+  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
+  /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.svg\"> <br>\n", stepm,YEARM,digitp,digit);
 */
-/*   fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.png\";replot;",digitp,optionfilefiname,digit); */
-  fprintf(ficgp,"\nset out \"%s%s.png\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
+/*   fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */
+  fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
 
   free_vector(xp,1,npar);
   free_matrix(doldm,1,nlstate,1,nlstate);
@@ -4156,19 +4206,19 @@ void varprob(char optionfilefiname[], do
   char fileresprobcor[FILENAMELENGTH];
   double ***varpij;
 
-  strcpy(fileresprob,"prob"); 
+  strcpy(fileresprob,"PROB_"); 
   strcat(fileresprob,fileres);
   if((ficresprob=fopen(fileresprob,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprob);
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
   }
-  strcpy(fileresprobcov,"probcov"); 
+  strcpy(fileresprobcov,"PROBCOV_"); 
   strcat(fileresprobcov,fileres);
   if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprobcov);
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
   }
-  strcpy(fileresprobcor,"probcor"); 
+  strcpy(fileresprobcor,"PROBCOR_"); 
   strcat(fileresprobcor,fileres);
   if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprobcor);
@@ -4211,10 +4261,9 @@ void varprob(char optionfilefiname[], do
   fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
   fprintf(fichtm,"\n");
 
-  fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of pairs of step probabilities (drawings)</a></h4></li>\n",optionfilehtmcov);
-  fprintf(fichtmcov,"\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n\
-  file %s<br>\n",optionfilehtmcov);
-  fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated\
+  fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.</li>\n",optionfilehtmcov);
+  fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
+  fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
 and drawn. It helps understanding how is the covariance between two incidences.\
  They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
   fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \
@@ -4235,23 +4284,23 @@ To be simple, these graphs help to under
     /*j1++;*/
       if  (cptcovn>0) {
 	fprintf(ficresprob, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresprob, "**********\n#\n");
 	fprintf(ficresprobcov, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresprobcov, "**********\n#\n");
 	
 	fprintf(ficgp, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficgp, "**********\n#\n");
 	
 	
 	fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
 	
 	fprintf(ficresprobcor, "\n#********** Variable ");    
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresprobcor, "**********\n#");    
       }
       
@@ -4264,7 +4313,8 @@ To be simple, these graphs help to under
 	if(nagesqr==1)
 	  cov[3]= age*age;
 	for (k=1; k<=cptcovn;k++) {
-	  cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4
+	  cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)];
+	  /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4
 							 * 1  1 1 1 1
 							 * 2  2 1 1 1
 							 * 3  1 2 1 1
@@ -4272,9 +4322,9 @@ To be simple, these graphs help to under
 	  /* nbcode[1][1]=0 nbcode[1][2]=1;*/
 	}
 	/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-	for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
+	for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
 	for (k=1; k<=cptcovprod;k++)
-	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
 	
     
 	for(theta=1; theta <=npar; theta++){
@@ -4422,17 +4472,18 @@ To be simple, these graphs help to under
 		  /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */
 		  if(first==1){
 		    first=0;
+ 		    fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n");
  		    fprintf(ficgp,"\nset parametric;unset label");
 		    fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
-		    fprintf(ficgp,"\nset ter png small size 320, 240");
+		    fprintf(ficgp,"\nset ter svg size 640, 480");
 		    fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
- :<a href=\"%s%d%1d%1d-%1d%1d.png\">\
-%s%d%1d%1d-%1d%1d.png</A>, ",k1,l1,k2,l2,\
-			    subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\
-			    subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
-		    fprintf(fichtmcov,"\n<br><img src=\"%s%d%1d%1d-%1d%1d.png\"> ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ :<a href=\"%s_%d%1d%1d-%1d%1d.svg\">\
+%s_%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\
+			    subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2,\
+			    subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+		    fprintf(fichtmcov,"\n<br><img src=\"%s_%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
 		    fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12);
-		    fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+		    fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
 		    fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
 		    fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
 		    fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\
@@ -4449,7 +4500,7 @@ To be simple, these graphs help to under
 		  }/* if first */
 		} /* age mod 5 */
 	      } /* end loop age */
-	      fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+	      fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
 	      first=1;
 	    } /*l12 */
 	  } /* k12 */
@@ -4471,7 +4522,7 @@ To be simple, these graphs help to under
 
 
 /******************* Printing html file ***********/
-void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
+void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
 		  int lastpass, int stepm, int weightopt, char model[],\
 		  int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
 		  int popforecast, int estepm ,\
@@ -4484,20 +4535,20 @@ void printinghtml(char fileres[], char t
 </ul>");
    fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \
  - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",
-	   jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"));
+	   jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));
    fprintf(fichtm,"\
  - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
-	   stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"));
+	   stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_"));
    fprintf(fichtm,"\
  - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
-	   subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
+	   subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));
    fprintf(fichtm,"\
  - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n",
-	   estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
+	   estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
    fprintf(fichtm,"\
  - Population projections by age and states: \
-   <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileres,"f"),subdirf2(fileres,"f"));
+   <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
 
 fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
 
@@ -4511,62 +4562,85 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++){ 
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
-	 printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);fflush(stdout);
+	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
+	 printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout);
        }
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }
+     /* aij, bij */
+     fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
+<img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
      /* Pij */
-     fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d_1.png\">%s%d_1.png</a><br> \
-<img src=\"%s%d_1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);     
+     fprintf(fichtm,"<br>\n- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
+<img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);     
      /* Quasi-incidences */
-     fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
- before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d_2.png\">%s%d_2.png</a><br> \
-<img src=\"%s%d_2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); 
-       /* Period (stable) prevalence in each health state */
-       for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \
-<img src=\"%s%d_%d.png\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
-       }
+     fprintf(fichtm,"<br>\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
+ before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too,\
+ incidence (rates) are the limit when h tends to zero of the ratio of the probability hPij \
+divided by h: hPij/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
+<img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); 
+     /* Survival functions (period) in state j */
+     for(cpt=1; cpt<=nlstate;cpt++){
+       fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1);
+     }
+     /* State specific survival functions (period) */
+     for(cpt=1; cpt<=nlstate;cpt++){
+       fprintf(fichtm,"<br>\n- Survival functions from state %d in any different live states and total.\
+ Or probability to survive in various states (1 to %d) being in state %d at different ages.\
+ <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1);
+     }
+     /* Period (stable) prevalence in each health state */
+     for(cpt=1; cpt<=nlstate;cpt++){
+       fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
+     }
      for(cpt=1; cpt<=nlstate;cpt++) {
-        fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
-<img src=\"%s%d%d.png\">",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
+       fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s%d%d.svg\">%s%d%d.svg</a> <br> \
+<img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
      }
    /* } /\* end i1 *\/ */
  }/* End k1 */
  fprintf(fichtm,"</ul>");
 
-
  fprintf(fichtm,"\
 \n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
- - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres);
+ - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \
+ - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.<br> \
+But because parameters are usually highly correlated (a higher incidence of disability \
+and a higher incidence of recovery can give very close observed transition) it might \
+be very useful to look not only at linear confidence intervals estimated from the \
+variances but at the covariance matrix. And instead of looking at the estimated coefficients \
+(parameters) of the logistic regression, it might be more meaningful to visualize the \
+covariance matrix of the one-step probabilities. \
+See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres);
 
- fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));
+ fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
+	 subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_"));
  fprintf(fichtm,"\
  - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov"));
+	 subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_"));
 
  fprintf(fichtm,"\
  - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));
+	 subdirf2(fileresu,"PROBCOR_"),subdirf2(fileresu,"PROBCOR_"));
  fprintf(fichtm,"\
  - Variances and covariances of health expectancies by age and <b>initial health status</b> (cov(e<sup>ij</sup>,e<sup>kl</sup>)(estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n</li>",
-	   estepm,subdirf2(fileres,"cve"),subdirf2(fileres,"cve"));
+	   estepm,subdirf2(fileresu,"CVE_"),subdirf2(fileresu,"CVE_"));
  fprintf(fichtm,"\
  - (a) Health expectancies by health status at initial age (e<sup>ij</sup>) and standard errors (in parentheses) (b) life expectancies and standard errors (e<sup>i.</sup>=e<sup>i1</sup>+e<sup>i2</sup>+...)(estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n</li>",
-	   estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde"));
+	   estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_"));
  fprintf(fichtm,"\
  - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
-	 estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
+	 estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_"));
  fprintf(fichtm,"\
  - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a> <br>\n",
-	 estepm, subdirf2(fileres,"t"),subdirf2(fileres,"t"));
+	 estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_"));
  fprintf(fichtm,"\
  - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
-	 subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
+	 subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
 
 /*  if(popforecast==1) fprintf(fichtm,"\n */
 /*  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */
@@ -4587,20 +4661,20 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++) 
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
+	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \
-prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.png <br>\
-<img src=\"%s%d_%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);  
+prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.svg <br>\
+<img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);  
      }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \
 health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
 true period expectancies (those weighted with period prevalences are also\
  drawn in addition to the population based expectancies computed using\
- observed and cahotic prevalences: %s%d.png<br>\
-<img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
+ observed and cahotic prevalences: %s_%d.svg<br>\
+<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);
    /* } /\* end i1 *\/ */
  }/* End k1 */
  fprintf(fichtm,"</ul>");
@@ -4608,11 +4682,12 @@ true period expectancies (those weighted
 }
 
 /******************* Gnuplot file **************/
-void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
 
   char dirfileres[132],optfileres[132];
   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
   int ng=0;
+  int vpopbased;
 /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
 /*     printf("Problem with file %s",optionfilegnuplot); */
 /*     fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */
@@ -4626,73 +4701,78 @@ void printinggnuplot(char fileres[], cha
   strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");
  /* 1eme*/
-  fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'vpl' files\n");
+  fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n");
   for (cpt=1; cpt<= nlstate ; cpt ++) {
     for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
-     fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1);
-     fprintf(ficgp,"\n#set out \"v%s%d_%d.png\" \n",optionfilefiname,cpt,k1);
+     fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
+     fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
      fprintf(ficgp,"set xlabel \"Age\" \n\
 set ylabel \"Probability\" \n\
-set ter png small size 320, 240\n\
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);
+set ter svg size 640, 480\n\
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
 
      for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else        fprintf(ficgp," %%*lf (%%*lf)");
      }
-     fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+     fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
      for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else fprintf(ficgp," %%*lf (%%*lf)");
      } 
-     fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); 
+     fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); 
      for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else fprintf(ficgp," %%*lf (%%*lf)");
      }  
-     fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
-   }
-  }
+     fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+     fprintf(ficgp,"\nset out \n");
+    } /* k1 */
+  } /* cpt */
   /*2 eme*/
   fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n");
   for (k1=1; k1<= m ; k1 ++) { 
-    fprintf(ficgp,"\nset out \"%s%d.png\" \n",subdirf2(optionfilefiname,"e"),k1);
-    fprintf(ficgp,"set ylabel \"Years\" \nset ter png small size 320, 240\nplot [%.f:%.f] ",ageminpar,fage);
-    
-    for (i=1; i<= nlstate+1 ; i ++) {
-      k=2*i;
-      fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
-      for (j=1; j<= nlstate+1 ; j ++) {
-	if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	else fprintf(ficgp," %%*lf (%%*lf)");
-      }   
-      if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
-      else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
-      fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
-      for (j=1; j<= nlstate+1 ; j ++) {
-	if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	else fprintf(ficgp," %%*lf (%%*lf)");
-      }   
-      fprintf(ficgp,"\" t\"\" w l lt 0,");
-      fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
-      for (j=1; j<= nlstate+1 ; j ++) {
-	if (j==i) fprintf(ficgp," %%lf (%%lf)");
-	else fprintf(ficgp," %%*lf (%%*lf)");
-      }   
-      if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
-      else fprintf(ficgp,"\" t\"\" w l lt 0,");
-    }
-  }
-  
+    fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
+    for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
+      if(vpopbased==0)
+	fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+      else
+	fprintf(ficgp,"\nreplot ");
+      for (i=1; i<= nlstate+1 ; i ++) {
+	k=2*i;
+	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
+	for (j=1; j<= nlstate+1 ; j ++) {
+	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	  else fprintf(ficgp," %%*lf (%%*lf)");
+	}   
+	if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
+	else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+	for (j=1; j<= nlstate+1 ; j ++) {
+	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	  else fprintf(ficgp," %%*lf (%%*lf)");
+	}   
+	fprintf(ficgp,"\" t\"\" w l lt 0,");
+	fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+	for (j=1; j<= nlstate+1 ; j ++) {
+	  if (j==i) fprintf(ficgp," %%lf (%%lf)");
+	  else fprintf(ficgp," %%*lf (%%*lf)");
+	}   
+	if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+	else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
+      } /* state */
+    } /* vpopbased */
+    fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
+  } /* k1 */
   /*3eme*/
   
   for (k1=1; k1<= m ; k1 ++) { 
     for (cpt=1; cpt<= nlstate ; cpt ++) {
       /*       k=2+nlstate*(2*cpt-2); */
       k=2+(nlstate+1)*(cpt-1);
-      fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1);
-      fprintf(ficgp,"set ter png small size 320, 240\n\
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileres,"e"),k1-1,k1-1,k,cpt);
+      fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
+      fprintf(ficgp,"set ter svg size 640, 480\n\
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
       /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
 	for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
 	fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
@@ -4702,39 +4782,98 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
 	
       */
       for (i=1; i< nlstate ; i ++) {
-	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+i,cpt,i+1);
+	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
 	/*	fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
 	
       } 
-      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+nlstate,cpt);
+      fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
     }
   }
   
+  /* Survival functions (period) from state i in state j by initial state i */
+  for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
+    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+      k=3;
+      fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'lij' files, cov=%d state=%d",k1, cpt);
+      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
+      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+set ter svg size 640, 480\n\
+unset log y\n\
+plot [%.f:%.f]  ", ageminpar, agemaxpar);
+      for (i=1; i<= nlstate ; i ++){
+	if(i==1)
+	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+	else
+	  fprintf(ficgp,", '' ");
+	l=(nlstate+ndeath)*(i-1)+1;
+	fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
+	for (j=2; j<= nlstate+ndeath ; j ++)
+	  fprintf(ficgp,"+$%d",k+l+j-1);
+	fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
+      } /* nlstate */
+      fprintf(ficgp,"\nset out\n");
+    } /* end cpt state*/ 
+  } /* end covariate */  
+
+  /* Survival functions (period) from state i in state j by final state j */
+  for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
+    for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
+      k=3;
+      fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
+      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
+      fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
+set ter svg size 640, 480\n\
+unset log y\n\
+plot [%.f:%.f]  ", ageminpar, agemaxpar);
+      for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+	if(j==1)
+	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+	else
+	  fprintf(ficgp,", '' ");
+	l=(nlstate+ndeath)*(cpt-1) +j;
+	fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
+	/* for (i=2; i<= nlstate+ndeath ; i ++) */
+	/*   fprintf(ficgp,"+$%d",k+l+i-1); */
+	fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
+      } /* nlstate */
+      fprintf(ficgp,", '' ");
+      fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
+      for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+	l=(nlstate+ndeath)*(cpt-1) +j;
+	if(j < nlstate)
+	  fprintf(ficgp,"$%d +",k+l);
+	else
+	  fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
+      }
+      fprintf(ficgp,"\nset out\n");
+    } /* end cpt state*/ 
+  } /* end covariate */  
+
   /* CV preval stable (period) */
   for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
       k=3;
       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);
-      fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1);
+      fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter png small size 320, 240\n\
+set ter svg size 640, 480\n\
 unset log y\n\
 plot [%.f:%.f]  ", ageminpar, agemaxpar);
       for (i=1; i<= nlstate ; i ++){
 	if(i==1)
-	  fprintf(ficgp,"\"%s\"",subdirf2(fileres,"pij"));
+	  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
 	else
 	  fprintf(ficgp,", '' ");
 	l=(nlstate+ndeath)*(i-1)+1;
 	fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
-	for (j=1; j<= (nlstate-1) ; j ++)
-	  fprintf(ficgp,"+$%d",k+l+j);
+	for (j=2; j<= nlstate ; j ++)
+	  fprintf(ficgp,"+$%d",k+l+j-1);
 	fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
       } /* nlstate */
-      fprintf(ficgp,"\n");
+      fprintf(ficgp,"\nset out\n");
     } /* end cpt state*/ 
   } /* end covariate */  
-  
+
   /* proba elementaires */
   fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
   for(i=1,jk=1; i <=nlstate; i++){
@@ -4753,7 +4892,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"##############\n#\n");
 
   /*goto avoid;*/
-  fprintf(ficgp,"\n##############\n#Graphics of of probabilities or incidences\n#############\n");
+  fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
   fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
   fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n");
   fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n");
@@ -4767,66 +4906,100 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   fprintf(ficgp,"#       +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n");
   fprintf(ficgp,"#       +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");
   fprintf(ficgp,"#\n");
-   for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
+   for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/
      fprintf(ficgp,"# ng=%d\n",ng);
      fprintf(ficgp,"#   jk=1 to 2^%d=%d\n",cptcoveff,m);
      for(jk=1; jk <=m; jk++) {
        fprintf(ficgp,"#    jk=%d\n",jk);
-       fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng); 
-       if (ng==2)
+       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);
+       fprintf(ficgp,"\nset ter svg size 640, 480 ");
+       if (ng==1){
+	 fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */
+	 fprintf(ficgp,"\nunset log y");
+       }else if (ng==2){
+	 fprintf(ficgp,"\nset ylabel \"Probability\"\n");
+	 fprintf(ficgp,"\nset log y");
+       }else if (ng==3){
 	 fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
-       else
-	 fprintf(ficgp,"\nset title \"Probability\"\n");
-       fprintf(ficgp,"\nset ter png small size 320, 240\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
+	 fprintf(ficgp,"\nset log y");
+       }else
+	 fprintf(ficgp,"\nunset title ");
+       fprintf(ficgp,"\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
        i=1;
        for(k2=1; k2<=nlstate; k2++) {
 	 k3=i;
 	 for(k=1; k<=(nlstate+ndeath); k++) {
 	   if (k != k2){
-	     if(ng==2)
+	     switch( ng) {
+	     case 1:
 	       if(nagesqr==0)
-		 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+		 fprintf(ficgp," p%d+p%d*x",i,i+1);
 	       else /* nagesqr =1 */
-		 fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
-	     else
+		 fprintf(ficgp," p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+	       break;
+	     case 2: /* ng=2 */
 	       if(nagesqr==0)
 		 fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
 	       else /* nagesqr =1 */
-		 fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+		   fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+	       break;
+	     case 3:
+	       if(nagesqr==0)
+		 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+	       else /* nagesqr =1 */
+		 fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
+	       break;
+	     }
 	     ij=1;/* To be checked else nbcode[0][0] wrong */
 	     for(j=3; j <=ncovmodel-nagesqr; j++) {
-	       if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */
-	       	 fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
-	       	 ij++;
+	       /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
+	       if(ij <=cptcovage) { /* Bug valgrind */
+		 if((j-2)==Tage[ij]) { /* Bug valgrind */
+		   fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+		   /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		   ij++;
+		 }
 	       }
 	       else
-		 fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+		 fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
 	     }
-	     fprintf(ficgp,")/(1");
+	     if(ng != 1){
+	       fprintf(ficgp,")/(1");
 	     
-	     for(k1=1; k1 <=nlstate; k1++){ 
-	       if(nagesqr==0)
-		 fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
-	       else /* nagesqr =1 */
-		 fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
-  
-	       ij=1;
-	       for(j=3; j <=ncovmodel-nagesqr; j++){
-		 if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
-		   fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
-		   ij++;
+	       for(k1=1; k1 <=nlstate; k1++){ 
+		 if(nagesqr==0)
+		   fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+		 else /* nagesqr =1 */
+		   fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
+		 
+		 ij=1;
+		 for(j=3; j <=ncovmodel-nagesqr; j++){
+		   if(ij <=cptcovage) { /* Bug valgrind */
+		     if((j-2)==Tage[ij]) { /* Bug valgrind */
+		       fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+		       /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+		       ij++;
+		     }
+		   }
+		   else
+		     fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
 		 }
-		 else
-		   fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+		 fprintf(ficgp,")");
 	       }
 	       fprintf(ficgp,")");
+	       if(ng ==2)
+		 fprintf(ficgp," t \"p%d%d\" ", k2,k);
+	       else /* ng= 3 */
+		 fprintf(ficgp," t \"i%d%d\" ", k2,k);
+	     }else{ /* end ng <> 1 */
+	       fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
 	     }
-	     fprintf(ficgp,") t \"p%d%d\" ", k2,k);
 	     if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
 	     i=i+ncovmodel;
 	   }
 	 } /* end k */
        } /* end k2 */
+       fprintf(ficgp,"\n set out\n");
      } /* end jk */
    } /* end ng */
  /* avoid: */
@@ -4894,8 +5067,8 @@ void prevforecast(char fileres[], double
   agelim=AGESUP;
   prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
  
-  strcpy(fileresf,"f"); 
-  strcat(fileresf,fileres);
+  strcpy(fileresf,"F_"); 
+  strcat(fileresf,fileresu);
   if((ficresf=fopen(fileresf,"w"))==NULL) {
     printf("Problem with forecast resultfile: %s\n", fileresf);
     fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
@@ -4944,7 +5117,7 @@ void prevforecast(char fileres[], double
       k=k+1;
       fprintf(ficresf,"\n#******");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
       fprintf(ficresf,"******\n");
       fprintf(ficresf,"# Covariate valuofcovar yearproj age");
@@ -4968,7 +5141,7 @@ void prevforecast(char fileres[], double
 	    if (h*hstepm/YEARM*stepm ==yearp) {
               fprintf(ficresf,"\n");
               for(j=1;j<=cptcoveff;j++) 
-                fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+                fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	      fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
 	    } 
 	    for(j=1; j<=nlstate+ndeath;j++) {
@@ -5018,8 +5191,8 @@ void populforecast(char fileres[], doubl
   prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
   
   
-  strcpy(filerespop,"pop"); 
-  strcat(filerespop,fileres);
+  strcpy(filerespop,"POP_"); 
+  strcat(filerespop,fileresu);
   if((ficrespop=fopen(filerespop,"w"))==NULL) {
     printf("Problem with forecast resultfile: %s\n", filerespop);
     fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop);
@@ -5066,7 +5239,7 @@ void populforecast(char fileres[], doubl
       k=k+1;
       fprintf(ficrespop,"\n#******");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
       fprintf(ficrespop,"******\n");
       fprintf(ficrespop,"# Age");
@@ -5372,7 +5545,7 @@ double gompertz_f(const gsl_vector *v, v
 #endif
 
 /******************* Printing html file ***********/
-void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \
+void printinghtmlmort(char fileresu[], char title[], char datafile[], int firstpass, \
 		  int lastpass, int stepm, int weightopt, char model[],\
 		  int imx,  double p[],double **matcov,double agemortsup){
   int i,k;
@@ -5381,7 +5554,7 @@ void printinghtmlmort(char fileres[], ch
   fprintf(fichtm,"  mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp);
   for (i=1;i<=2;i++) 
     fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
-  fprintf(fichtm,"<br><br><img src=\"graphmort.png\">");
+  fprintf(fichtm,"<br><br><img src=\"graphmort.svg\">");
   fprintf(fichtm,"</ul>");
 
 fprintf(fichtm,"<ul><li><h4>Life table</h4>\n <br>");
@@ -5396,7 +5569,7 @@ fprintf(fichtm,"<ul><li><h4>Life table</
 }
 
 /******************* Gnuplot file **************/
-void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+void printinggnuplotmort(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
 
   char dirfileres[132],optfileres[132];
 
@@ -5410,9 +5583,9 @@ void printinggnuplotmort(char fileres[],
 
   strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");
-  fprintf(ficgp,"set out \"graphmort.png\"\n "); 
+  fprintf(ficgp,"set out \"graphmort.svg\"\n "); 
   fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n "); 
-  fprintf(ficgp, "set ter png small size 320, 240\n set log y\n"); 
+  fprintf(ficgp, "set ter svg size 640, 480\n set log y\n"); 
   /* fprintf(ficgp, "set size 0.65,0.65\n"); */
   fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp);
 
@@ -5433,8 +5606,8 @@ int readdata(char datafile[], int firsto
 
 
   if((fic=fopen(datafile,"r"))==NULL)    {
-    printf("Problem while opening datafile: %s\n", datafile);return 1;
-    fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);return 1;
+    printf("Problem while opening datafile: %s\n", datafile);fflush(stdout);
+    fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);fflush(ficlog);return 1;
   }
 
   i=1;
@@ -5723,9 +5896,9 @@ int decodemodel ( char model[], int last
       /* 	k=1 Tvar[1]=2 (from V2) */
       /* 	k=5 Tvar[5] */
       /* for (k=1; k<=cptcovn;k++) { */
-      /* 	cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */
+      /* 	cov[2+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
       /* 	} */
-      /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; */
+      /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k])]]*cov[2]; */
       /*
        * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */
       for(k=cptcovt; k>=1;k--) /**< Number of covariates */
@@ -6109,8 +6282,8 @@ int prevalence_limit(double *p, double *
   double ftolpl = 1.e-10;
   double age, agebase, agelim;
 
-    strcpy(filerespl,"pl");
-    strcat(filerespl,fileres);
+    strcpy(filerespl,"PL_");
+    strcat(filerespl,fileresu);
     if((ficrespl=fopen(filerespl,"w"))==NULL) {
       printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
       fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
@@ -6136,14 +6309,14 @@ int prevalence_limit(double *p, double *
       //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
 	k=k+1;
 	/* to clean */
-	//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]);
-	fprintf(ficrespl,"\n#******");
-	printf("\n#******");
-	fprintf(ficlog,"\n#******");
+	//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+	fprintf(ficrespl,"#******");
+	printf("#******");
+	fprintf(ficlog,"#******");
 	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-	  printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-	  fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	  printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	  fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	fprintf(ficrespl,"******\n");
 	printf("******\n");
@@ -6151,7 +6324,7 @@ int prevalence_limit(double *p, double *
 
 	fprintf(ficrespl,"#Age ");
 	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
 	fprintf(ficrespl,"\n");
@@ -6161,7 +6334,7 @@ int prevalence_limit(double *p, double *
 	  prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
 	  fprintf(ficrespl,"%.0f ",age );
 	  for(j=1;j<=cptcoveff;j++)
-	    fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	    fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	  for(i=1; i<=nlstate;i++)
 	    fprintf(ficrespl," %.5f", prlim[i][i]);
 	  fprintf(ficrespl,"\n");
@@ -6183,7 +6356,7 @@ int hPijx(double *p, int bage, int fage)
   double agedeb;
   double ***p3mat;
 
-    strcpy(filerespij,"pij");  strcat(filerespij,fileres);
+    strcpy(filerespij,"PIJ_");  strcat(filerespij,fileresu);
     if((ficrespij=fopen(filerespij,"w"))==NULL) {
       printf("Problem with Pij resultfile: %s\n", filerespij); return 1;
       fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij); return 1;
@@ -6208,7 +6381,7 @@ int hPijx(double *p, int bage, int fage)
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
       fprintf(ficrespij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) 
-	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       fprintf(ficrespij,"******\n");
       
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
@@ -6260,6 +6433,7 @@ int main(int argc, char *argv[])
 
   int jj, ll, li, lj, lk;
   int numlinepar=0; /* Current linenumber of parameter file */
+  int num_filled;
   int itimes;
   int NDIM=2;
   int vpopbased=0;
@@ -6269,7 +6443,8 @@ int main(int argc, char *argv[])
   /* FILE *ficgp;*/ /*Gnuplot File */
   struct stat info;
   double agedeb=0.;
-  double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
+
+  double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
 
   double fret;
   double dum=0.; /* Dummy variable */
@@ -6277,13 +6452,17 @@ int main(int argc, char *argv[])
   double ***mobaverage;
 
   char line[MAXLINE];
-  char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
+  char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
+
+  char model[MAXLINE], modeltemp[MAXLINE];
   char pathr[MAXLINE], pathimach[MAXLINE]; 
   char *tok, *val; /* pathtot */
   int firstobs=1, lastobs=10;
-  int c,  h , cpt;
+  int c,  h , cpt, c2;
   int jl=0;
   int i1, j1, jk, stepsize=0;
+  int count=0;
+
   int *tab; 
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
   int mobilav=0,popforecast=0;
@@ -6355,7 +6534,7 @@ int main(int argc, char *argv[])
   getcwd(pathcd, size);
 #endif
   syscompilerinfo(0);
-  printf("\n%s\n%s",version,fullversion);
+  printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion);
   if(argc <=1){
     printf("\nEnter the parameter file name: ");
     fgets(pathr,FILENAMELENGTH,stdin);
@@ -6419,7 +6598,7 @@ int main(int argc, char *argv[])
     goto end;
   }
   fprintf(ficlog,"Log filename:%s\n",filelog);
-  fprintf(ficlog,"\n%s\n%s",version,fullversion);
+  fprintf(ficlog,"Version %s %s",version,fullversion);
   fprintf(ficlog,"\nEnter the parameter file name: \n");
   fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
  path=%s \n\
@@ -6427,7 +6606,7 @@ int main(int argc, char *argv[])
  optionfilext=%s\n\
  optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
 
-  syscompilerinfo(0);
+  syscompilerinfo(1);
 
   printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Local time (at start): %s",strstart);
@@ -6438,7 +6617,9 @@ int main(int argc, char *argv[])
   /* */
   strcpy(fileres,"r");
   strcat(fileres, optionfilefiname);
+  strcat(fileresu, optionfilefiname); /* Without r in front */
   strcat(fileres,".txt");    /* Other files have txt extension */
+  strcat(fileresu,".txt");    /* Other files have txt extension */
 
   /* Main ---------arguments file --------*/
 
@@ -6453,7 +6634,7 @@ int main(int argc, char *argv[])
 
 
   strcpy(filereso,"o");
-  strcat(filereso,fileres);
+  strcat(filereso,fileresu);
   if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */
     printf("Problem with Output resultfile: %s\n", filereso);
     fprintf(ficlog,"Problem with Output resultfile: %s\n", filereso);
@@ -6463,23 +6644,84 @@ int main(int argc, char *argv[])
 
   /* Reads comments: lines beginning with '#' */
   numlinepar=0;
-  while((c=getc(ficpar))=='#' && c!= EOF){
-    ungetc(c,ficpar);
-    fgets(line, MAXLINE, ficpar);
+
+    /* First parameter line */
+  while(fgets(line, MAXLINE, ficpar)) {
+    /* If line starts with a # it is a comment */
+    if (line[0] == '#') {
+      numlinepar++;
+      fputs(line,stdout);
+      fputs(line,ficparo);
+      fputs(line,ficlog);
+      continue;
+    }else
+      break;
+  }
+  if((num_filled=sscanf(line,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", \
+			title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){
+    if (num_filled != 5) {
+      printf("Should be 5 parameters\n");
+    }
     numlinepar++;
-    fputs(line,stdout);
-    fputs(line,ficparo);
-    fputs(line,ficlog);
+    printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass);
+  }
+  /* Second parameter line */
+  while(fgets(line, MAXLINE, ficpar)) {
+    /* If line starts with a # it is a comment */
+    if (line[0] == '#') {
+      numlinepar++;
+      fputs(line,stdout);
+      fputs(line,ficparo);
+      fputs(line,ficlog);
+      continue;
+    }else
+      break;
+  }
+  if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
+			&ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
+    if (num_filled != 8) {
+      printf("Not 8\n");
+    }
+    printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
   }
-  ungetc(c,ficpar);
 
-  fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
-  numlinepar++;
-  printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
+  /* Third parameter line */
+  while(fgets(line, MAXLINE, ficpar)) {
+    /* If line starts with a # it is a comment */
+    if (line[0] == '#') {
+      numlinepar++;
+      fputs(line,stdout);
+      fputs(line,ficparo);
+      fputs(line,ficlog);
+      continue;
+    }else
+      break;
+  }
+  if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
+    if (num_filled == 0)
+            model[0]='\0';
+    else if (num_filled != 1){
+      printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
+      fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
+      model[0]='\0';
+      goto end;
+    }
+    else{
+      if (model[0]=='+'){
+	for(i=1; i<=strlen(model);i++)
+	  modeltemp[i-1]=model[i];
+	strcpy(model,modeltemp); 
+      }
+    }
+    /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */
+  }
+  /* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */
+  /* numlinepar=numlinepar+3; /\* In general *\/ */
+  /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
   if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */
     model[strlen(model)-1]='\0';
-  fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
-  fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+  fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+  fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
   fflush(ficlog);
   /* if(model[0]=='#'|| model[0]== '\0'){ */
   if(model[0]=='#'){
@@ -6495,6 +6737,10 @@ int main(int argc, char *argv[])
     ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);
     numlinepar++;
+    if(line[1]=='q'){ /* This #q will quit imach (the answer is q) */
+      z[0]=line[1];
+    }
+    /* printf("****line [1] = %c \n",line[1]); */
     fputs(line, stdout);
     //puts(line);
     fputs(line,ficparo);
@@ -6561,7 +6807,7 @@ int main(int argc, char *argv[])
 	if(jj==i) continue;
 	j++;
 	fscanf(ficpar,"%1d%1d",&i1,&j1);
-	if ((i1 != i) && (j1 != j)){
+	if ((i1 != i) || (j1 != jj)){
 	  printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
 It might be a problem of design; if ncovcol and the model are correct\n \
 run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);
@@ -6569,8 +6815,8 @@ run imach with mle=-1 to get a correct t
 	}
 	fprintf(ficparo,"%1d%1d",i1,j1);
 	if(mle==1)
-	  printf("%1d%1d",i,j);
-	fprintf(ficlog,"%1d%1d",i,j);
+	  printf("%1d%1d",i,jj);
+	fprintf(ficlog,"%1d%1d",i,jj);
 	for(k=1; k<=ncovmodel;k++){
 	  fscanf(ficpar," %lf",&param[i][j][k]);
 	  if(mle==1){
@@ -6651,12 +6897,22 @@ run imach with mle=-1 to get a correct t
     for(i=1; i <=npar; i++)
       for(j=1; j <=npar; j++) matcov[i][j]=0.;
       
+    /* Scans npar lines */
     for(i=1; i <=npar; i++){
-      fscanf(ficpar,"%s",str);
+      count=fscanf(ficpar,"%1d%1d%1d",&i1,&j1,&jk);
+      if(count != 3){
+	printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
+This is probably because your covariance matrix doesn't \n  contain exactly %d lines corresponding to your model line '1+age+%s'.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
+	fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
+This is probably because your covariance matrix doesn't \n  contain exactly %d lines corresponding to your model line '1+age+%s'.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
+	exit(1);
+      }else
       if(mle==1)
-	printf("%s",str);
-      fprintf(ficlog,"%s",str);
-      fprintf(ficparo,"%s",str);
+	printf("%1d%1d%1d",i1,j1,jk);
+      fprintf(ficlog,"%1d%1d%1d",i1,j1,jk);
+      fprintf(ficparo,"%1d%1d%1d",i1,j1,jk);
       for(j=1; j <=i; j++){
 	fscanf(ficpar," %le",&matcov[i][j]);
 	if(mle==1){
@@ -6672,6 +6928,7 @@ run imach with mle=-1 to get a correct t
       fprintf(ficlog,"\n");
       fprintf(ficparo,"\n");
     }
+    /* End of read covariance matrix npar lines */
     for(i=1; i <=npar; i++)
       for(j=i+1;j<=npar;j++)
 	matcov[i][j]=matcov[j][i];
@@ -6688,8 +6945,8 @@ run imach with mle=-1 to get a correct t
     strcat(rfileres,".");    /* */
     strcat(rfileres,optionfilext);    /* Other files have txt extension */
     if((ficres =fopen(rfileres,"w"))==NULL) {
-      printf("Problem writing new parameter file: %s\n", fileres);goto end;
-      fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
+      printf("Problem writing new parameter file: %s\n", rfileres);goto end;
+      fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end;
     }
     fprintf(ficres,"#%s\n",version);
   }    /* End of mle != -3 */
@@ -6799,8 +7056,8 @@ run imach with mle=-1 to get a correct t
      V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
   /* 1 to ncodemax[j] is the maximum value of this jth covariate */
 
-  codtab=imatrix(1,100,1,10); /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
-  /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtab[100][10]);*/
+  /*  codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
+  /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
   /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
   h=0;
 
@@ -6810,13 +7067,6 @@ run imach with mle=-1 to get a correct t
  
   m=pow(2,cptcoveff);
  
-  for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */
-    for(i=1; i <=pow(2,cptcoveff-k);i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ 
-      for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate ncodemax=2*/
-	for(cpt=1; cpt <=pow(2,k-1); cpt++){  /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ 
-	  h++;
-	  if (h>m) 
-	    h=1;
 	  /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
 	   * For k=4 covariates, h goes from 1 to 2**k
 	   * codtabm(h,k)=  1 & (h-1) >> (k-1) ;
@@ -6830,32 +7080,47 @@ run imach with mle=-1 to get a correct t
 	   *     6     2     1     2     1
 	   *     7 i=4 1     2     2     1
 	   *     8     2     2     2     1
-	   *     9 i=5 1 i=3 1 i=2 1     1
-	   *    10     2     1     1     1
-	   *    11 i=6 1     2     1     1
-	   *    12     2     2     1     1
-	   *    13 i=7 1 i=4 1     2     1    
-	   *    14     2     1     2     1
-	   *    15 i=8 1     2     2     1
-	   *    16     2     2     2     1
+	   *     9 i=5 1 i=3 1 i=2 1     2
+	   *    10     2     1     1     2
+	   *    11 i=6 1     2     1     2
+	   *    12     2     2     1     2
+	   *    13 i=7 1 i=4 1     2     2    
+	   *    14     2     1     2     2
+	   *    15 i=8 1     2     2     2
+	   *    16     2     2     2     2
 	   */
-	  codtab[h][k]=j;
-	  /* codtab[12][3]=1; */
-	  /*codtab[h][Tvar[k]]=j;*/
-	  printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]);
-	} 
-      }
-    }
-  } 
+  for(h=1; h <=100 ;h++){ 
+    /* printf("h=%2d ", h); */
+     /* for(k=1; k <=10; k++){ */
+       /* printf("k=%d %d ",k,codtabm(h,k)); */
+     /*   codtab[h][k]=codtabm(h,k); */
+     /* } */
+     /* printf("\n"); */
+  }
+  /* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */
+  /*   for(i=1; i <=pow(2,cptcoveff-k);i++){ /\* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 *\/  */
+  /*     for(j=1; j <= ncodemax[k]; j++){ /\* For each modality of this covariate ncodemax=2*\/ */
+  /* 	for(cpt=1; cpt <=pow(2,k-1); cpt++){  /\* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 *\/  */
+  /* 	  h++; */
+  /* 	  if (h>m)  */
+  /* 	    h=1; */
+  /* 	  codtab[h][k]=j; */
+  /* 	  /\* codtab[12][3]=1; *\/ */
+  /* 	  /\*codtab[h][Tvar[k]]=j;*\/ */
+  /* 	  /\* printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); *\/ */
+  /* 	}  */
+  /*     } */
+  /*   } */
+  /* }  */
   /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); 
      codtab[1][2]=1;codtab[2][2]=2; */
-  /* for(i=1; i <=m ;i++){ 
-     for(k=1; k <=cptcovn; k++){
-       printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
-     }
-     printf("\n");
-     }
-     scanf("%d",i);*/
+  /* for(i=1; i <=m ;i++){  */
+  /*    for(k=1; k <=cptcovn; k++){ */
+  /*      printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); */
+  /*    } */
+  /*    printf("\n"); */
+  /* } */
+  /*   scanf("%d",i);*/
 
  free_ivector(Ndum,-1,NCOVMAX);
 
@@ -6864,7 +7129,7 @@ run imach with mle=-1 to get a correct t
   /* Initialisation of ----------- gnuplot -------------*/
   strcpy(optionfilegnuplot,optionfilefiname);
   if(mle==-3)
-    strcat(optionfilegnuplot,"-mort");
+    strcat(optionfilegnuplot,"-MORT_");
   strcat(optionfilegnuplot,".gp");
 
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
@@ -6883,7 +7148,7 @@ run imach with mle=-1 to get a correct t
 
   strcpy(optionfilehtm,optionfilefiname); /* Main html file */
   if(mle==-3)
-    strcat(optionfilehtm,"-mort");
+    strcat(optionfilehtm,"-MORT_");
   strcat(optionfilehtm,".htm");
   if((fichtm=fopen(optionfilehtm,"w"))==NULL)    {
     printf("Problem with %s \n",optionfilehtm);
@@ -6999,8 +7264,8 @@ Interval (in months) between two waves:
 #else
     printf("Powell\n");  fprintf(ficlog,"Powell\n");
 #endif
-    strcpy(filerespow,"pow-mort"); 
-    strcat(filerespow,fileres);
+    strcpy(filerespow,"POW-MORT_"); 
+    strcat(filerespow,fileresu);
     if((ficrespow=fopen(filerespow,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", filerespow);
       fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
@@ -7111,9 +7376,10 @@ Interval (in months) between two waves:
     }
     
     printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);
-    for (i=1;i<=NDIM;i++) 
+    for (i=1;i<=NDIM;i++) {
       printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
-
+      fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
+    }
     lsurv=vector(1,AGESUP);
     lpop=vector(1,AGESUP);
     tpop=vector(1,AGESUP);
@@ -7145,9 +7411,16 @@ Interval (in months) between two waves:
     
     
     replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
-    printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
-    
-    printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
+    if(ageminpar == AGEOVERFLOW ||agemaxpar == AGEOVERFLOW){
+      	printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
+This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
+      	fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
+This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
+    }else
+      printinggnuplotmort(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+    printinghtmlmort(fileresu,title,datafile, firstpass, lastpass, \
 		     stepm, weightopt,\
 		     model,imx,p,matcov,agemortsup);
     
@@ -7210,6 +7483,24 @@ Interval (in months) between two waves:
       ftolhess=ftol; /* Usually correct */
       hesscov(matcov, p, npar, delti, ftolhess, func);
     }
+    printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
+    fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n  It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
+    for(i=1,jk=1; i <=nlstate; i++){
+      for(k=1; k <=(nlstate+ndeath); k++){
+	if (k != i) {
+	  printf("%d%d ",i,k);
+	  fprintf(ficlog,"%d%d ",i,k);
+	  for(j=1; j <=ncovmodel; j++){
+	    printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
+	    fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
+	    jk++; 
+	  }
+	  printf("\n");
+	  fprintf(ficlog,"\n");
+	}
+      }
+    }
+
     fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
     printf("# Scales (for hessian or gradient estimation)\n");
     fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
@@ -7366,6 +7657,7 @@ Interval (in months) between two waves:
     dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
     
     fscanf(ficpar,"pop_based=%d\n",&popbased);
+    fprintf(ficlog,"pop_based=%d\n",popbased);
     fprintf(ficparo,"pop_based=%d\n",popbased);   
     fprintf(ficres,"pop_based=%d\n",popbased);   
     
@@ -7390,9 +7682,17 @@ Interval (in months) between two waves:
     /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */
     
     replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
-    printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+    if(ageminpar == AGEOVERFLOW ||agemaxpar == -AGEOVERFLOW){
+      	printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
+This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
+      	fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
+This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
+    }else
+      printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
     
-    printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\
+    printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
 		 model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
 		 jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
       
@@ -7444,7 +7744,7 @@ Interval (in months) between two waves:
     /*if((stepm == 1) && (strcmp(model,".")==0)){*/
     if(prevfcast==1){
       /*    if(stepm ==1){*/
-      prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
+      prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
       /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
       /*      }  */
       /*      else{ */
@@ -7474,8 +7774,8 @@ Interval (in months) between two waves:
 
     /*---------- Health expectancies, no variances ------------*/
 
-    strcpy(filerese,"e");
-    strcat(filerese,fileres);
+    strcpy(filerese,"E_");
+    strcat(filerese,fileresu);
     if((ficreseij=fopen(filerese,"w"))==NULL) {
       printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
       fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
@@ -7488,7 +7788,7 @@ Interval (in months) between two waves:
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
 	fprintf(ficreseij,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	fprintf(ficreseij,"******\n");
 
@@ -7505,8 +7805,8 @@ Interval (in months) between two waves:
     /*---------- Health expectancies and variances ------------*/
 
 
-    strcpy(filerest,"t");
-    strcat(filerest,fileres);
+    strcpy(filerest,"T_");
+    strcat(filerest,fileresu);
     if((ficrest=fopen(filerest,"w"))==NULL) {
       printf("Problem with total LE resultfile: %s\n", filerest);goto end;
       fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
@@ -7515,8 +7815,8 @@ Interval (in months) between two waves:
     fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); 
 
 
-    strcpy(fileresstde,"stde");
-    strcat(fileresstde,fileres);
+    strcpy(fileresstde,"STDE_");
+    strcat(fileresstde,fileresu);
     if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {
       printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
       fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
@@ -7524,8 +7824,8 @@ Interval (in months) between two waves:
     printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
     fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
 
-    strcpy(filerescve,"cve");
-    strcat(filerescve,fileres);
+    strcpy(filerescve,"CVE_");
+    strcat(filerescve,fileresu);
     if((ficrescveij=fopen(filerescve,"w"))==NULL) {
       printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
       fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
@@ -7533,8 +7833,8 @@ Interval (in months) between two waves:
     printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
     fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
 
-    strcpy(fileresv,"v");
-    strcat(fileresv,fileres);
+    strcpy(fileresv,"V_");
+    strcat(fileresv,fileresu);
     if((ficresvij=fopen(fileresv,"w"))==NULL) {
       printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
       fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
@@ -7548,21 +7848,21 @@ Interval (in months) between two waves:
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
     	fprintf(ficrest,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) 
-	  fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	fprintf(ficrest,"******\n");
 
 	fprintf(ficresstdeij,"\n#****** ");
 	fprintf(ficrescveij,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-	  fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	  fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	fprintf(ficresstdeij,"******\n");
 	fprintf(ficrescveij,"******\n");
 
 	fprintf(ficresvij,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) 
-	  fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	fprintf(ficresvij,"******\n");
 
 	eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
@@ -7577,7 +7877,7 @@ Interval (in months) between two waves:
 
 
 	for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-	  oldm=oldms;savm=savms; /* Segmentation fault */
+	  oldm=oldms;savm=savms; /* ZZ Segmentation fault */
 	  cptcod= 0; /* To be deleted */
 	  varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
 	  fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
@@ -7585,13 +7885,13 @@ Interval (in months) between two waves:
 	    fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
 	  else
 	    fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
-	  fprintf(ficrest,"# Age e.. (std) ");
+	  fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
 	  for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
 	  fprintf(ficrest,"\n");
-
+	  /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
 	  epj=vector(1,nlstate+1);
 	  for(age=bage; age <=fage ;age++){
-	    prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+	    prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); /*ZZ Is it the correct prevalim */
 	    if (vpopbased==1) {
 	      if(mobilav ==0){
 		for(i=1; i<=nlstate;i++)
@@ -7602,14 +7902,17 @@ Interval (in months) between two waves:
 	      }
 	    }
 	
-	    fprintf(ficrest," %4.0f",age);
+	    fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
+	    /* printf(" age %4.0f ",age); */
 	    for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
 	      for(i=1, epj[j]=0.;i <=nlstate;i++) {
 		epj[j] += prlim[i][i]*eij[i][j][(int)age];
-		/*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+		/*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+	        /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
 	      }
 	      epj[nlstate+1] +=epj[j];
 	    }
+	    /* printf(" age %4.0f \n",age); */
 
 	    for(i=1, vepp=0.;i <=nlstate;i++)
 	      for(j=1;j <=nlstate;j++)
@@ -7641,8 +7944,8 @@ Interval (in months) between two waves:
   
     /*------- Variance of period (stable) prevalence------*/   
 
-    strcpy(fileresvpl,"vpl");
-    strcat(fileresvpl,fileres);
+    strcpy(fileresvpl,"VPL_");
+    strcat(fileresvpl,fileresu);
     if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
       printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);
       exit(0);
@@ -7655,7 +7958,7 @@ Interval (in months) between two waves:
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
     	fprintf(ficresvpl,"\n#****** ");
 	for(j=1;j<=cptcoveff;j++) 
-	  fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	fprintf(ficresvpl,"******\n");
       
 	varpl=matrix(1,nlstate,(int) bage, (int) fage);
@@ -7692,7 +7995,7 @@ Interval (in months) between two waves:
     free_ivector(Tage,1,NCOVMAX);
 
     free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
-    free_imatrix(codtab,1,100,1,10);
+    /* free_imatrix(codtab,1,100,1,10); */
   fflush(fichtm);
   fflush(ficgp);
   
@@ -7793,7 +8096,7 @@ Interval (in months) between two waves:
   }
   end:
   while (z[0] != 'q') {
-    printf("\nType  q for exiting: ");
+    printf("\nType  q for exiting: "); fflush(stdout);
     scanf("%s",z);
   }
 }