--- imach/src/imach.c	2015/09/09 16:53:55	1.200
+++ imach/src/imach.c	2015/09/15 17:34:58	1.201
@@ -1,6 +1,12 @@
-/* $Id: imach.c,v 1.200 2015/09/09 16:53:55 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
 
@@ -725,12 +731,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.200 2015/09/09 16:53:55 brouard Exp $ */
+/* $Id: imach.c,v 1.201 2015/09/15 17:34:58 brouard Exp $ */
 /* $State: Exp $ */
 #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.200 $ $Date: 2015/09/09 16:53:55 $"; 
+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  */
@@ -796,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];
@@ -2570,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);
@@ -2587,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;
@@ -2620,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);
@@ -2947,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);
@@ -3816,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);
@@ -3830,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 */
@@ -4070,7 +4076,7 @@ void varevsij(char optionfilefiname[], d
   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,"\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); */
@@ -4078,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.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",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.svg\";replot;",digitp,optionfilefiname,digit); */
-  fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),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);
@@ -4200,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);
@@ -4471,13 +4477,13 @@ To be simple, these graphs help to under
 		    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 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.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);
+ :<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.svg\"",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",\
@@ -4494,7 +4500,7 @@ To be simple, these graphs help to under
 		  }/* if first */
 		} /* age mod 5 */
 	      } /* end loop age */
-	      fprintf(ficgp,"\nset out;\nset out \"%s%d%1d%1d-%1d%1d.svg\";replot;set out;",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 */
@@ -4516,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 ,\
@@ -4529,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>");
 
@@ -4561,21 +4567,37 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
        }
        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.svg\">%s%d_1.svg</a><br> \
-<img src=\"%s%d_1.svg\">",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.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); 
-       /* 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.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);
-       }
+     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.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);
+       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 */
@@ -4594,31 +4616,31 @@ covariance matrix of the one-step probab
 See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres);
 
  fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
-	 subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));
+	 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 */
@@ -4645,14 +4667,14 @@ See page 'Matrix of variance-covariance
      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.svg <br>\
-<img src=\"%s%d_%d.svg\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);  
+<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.svg<br>\
-<img src=\"%s%d.svg\">",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>");
@@ -4660,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); */
@@ -4678,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.svg\" \n",subdirf2(optionfilefiname,"v"),cpt,k1);
-     fprintf(ficgp,"\n#set out \"v%s%d_%d.svg\" \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 svg size 640, 480\n\
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);
+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.svg\" \n",subdirf2(optionfilefiname,"e"),k1);
-    fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\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.svg\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1);
+      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(fileres,"e"),k1-1,k1-1,k,cpt);
+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);
@@ -4754,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.svg\" \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 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++){
@@ -4819,32 +4906,50 @@ 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.svg\" \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,"\nunset title \n");
-       fprintf(ficgp,"\nset ter svg size 640, 480\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++) {
 	       /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
@@ -4858,34 +4963,43 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 	       else
 		 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(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++;
+	       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]][codtabm(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: */
@@ -4953,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);
@@ -5077,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);
@@ -5431,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;
@@ -5455,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];
 
@@ -6168,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;
@@ -6242,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;
@@ -6503,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 --------*/
 
@@ -6518,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);
@@ -6581,8 +6697,10 @@ int main(int argc, char *argv[])
     }else
       break;
   }
-  if((num_filled=sscanf(line,"model=1+age%[^.\n]\n", model)) !=EOF){
-    if (num_filled != 1) {
+  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';
@@ -6592,8 +6710,8 @@ int main(int argc, char *argv[])
       if (model[0]=='+'){
 	for(i=1; i<=strlen(model);i++)
 	  modeltemp[i-1]=model[i];
+	strcpy(model,modeltemp); 
       }
-      strcpy(model,modeltemp); 
     }
     /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */
   }
@@ -6827,8 +6945,8 @@ Please run with mle=-1 to get a correct
     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 */
@@ -7011,7 +7129,7 @@ Please run with mle=-1 to get a correct
   /* Initialisation of ----------- gnuplot -------------*/
   strcpy(optionfilegnuplot,optionfilefiname);
   if(mle==-3)
-    strcat(optionfilegnuplot,"-mort");
+    strcat(optionfilegnuplot,"-MORT_");
   strcat(optionfilegnuplot,".gp");
 
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
@@ -7030,7 +7148,7 @@ Please run with mle=-1 to get a correct
 
   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);
@@ -7146,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);
@@ -7301,8 +7419,8 @@ Please run with mle=-1 to get a correct
 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(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
-    printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
+      printinggnuplotmort(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+    printinghtmlmort(fileresu,title,datafile, firstpass, lastpass, \
 		     stepm, weightopt,\
 		     model,imx,p,matcov,agemortsup);
     
@@ -7572,9 +7690,9 @@ Please run with mle=-1 to get a correct
 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(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+      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);
       
@@ -7626,7 +7744,7 @@ Please run with mle=-1 to get a correct
     /*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{ */
@@ -7656,8 +7774,8 @@ Please run with mle=-1 to get a correct
 
     /*---------- 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);
@@ -7687,8 +7805,8 @@ Please run with mle=-1 to get a correct
     /*---------- 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;
@@ -7697,8 +7815,8 @@ Please run with mle=-1 to get a correct
     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);
@@ -7706,8 +7824,8 @@ Please run with mle=-1 to get a correct
     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);
@@ -7715,8 +7833,8 @@ Please run with mle=-1 to get a correct
     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);
@@ -7767,7 +7885,7 @@ Please run with mle=-1 to get a correct
 	    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"); */
@@ -7784,7 +7902,7 @@ Please run with mle=-1 to get a correct
 	      }
 	    }
 	
-	    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++) {
@@ -7826,8 +7944,8 @@ Please run with mle=-1 to get a correct
   
     /*------- 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);