#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#include <assert.h>
#include <X11/Intrinsic.h>
#include <X11/StringDefs.h>
#include <X11/Shell.h>
#include <X11/Xlib.h>
#include <X11/Xaw/Cardinals.h>
#include <X11/Xaw/Command.h>
#include <X11/Xaw/Label.h>
#include <X11/Xaw/Box.h>
#include <X11/Xaw/List.h>
#include <X11/Xaw/Scrollbar.h>
#include <X11/Xaw/Dialog.h>
#include <X11/Xaw/Simple.h>
#include <X11/Xaw/Form.h>
#include <X11/Xaw/AsciiText.h>

#define Xmax 		10.		// spatial interval: from 0 to Xmax (mm)
#define N_U			1000	// Number of discrete points of a grid for U
#define N_V			100		// Number of discrete points of a grid for V
#define x_dermis	2.		// The width of the dermis region (mm)
#define epsilon		0.03	//	mm
#define m0			0.1
#define D_u			0.01	//	mm^2/h
#define mu_m		0.05	//	1/h
#define GAMMA		0.1
#define V_star		10.
#define k_d			0.01	//	1/h
#define V_0			1.
#define TIME_final	100.	//	h
typedef double TYPE;
typedef TYPE *ARRAY;
ARRAY X_u,X_v,U_n,U_n_1,gradU_n,gradU_n_1,V_n,V_n_1,AU,BU,CU,DU;

typedef int TYPE1;
typedef TYPE1 *INDEX;
INDEX KILL,DIVIDE;

struct cells {
double x;
double y;   // this coordinate is used only for visualization
double rmax;
double pmax;
double emax;
double r_i;
double e_i;
double delta_m;
double nabla_m;
};
struct cells *CELL_n,*CELL_n_1;

struct TEMP {
double x;
double y;   // this coordinate is used only for visualization
double rmax;
double pmax;
double emax;
};
struct TEMP *temp;

void ModelParameters(),SolvePDE(),CellDensity(),ComputeGradient(double U[],double gU[]),MoveCells(),
       InitialRandomSeed(),PrintMovieData(),Create_Arrays(),Average_and_print_depth();
double Tridiag(double a[],double b[],double c[],double f[],double d[]),
       ReceptorLevel(double U[], double x_i, double rmax_i),Average_r_max(int J),
       Average_p_max(int J),Average_e_max(int J),Average_Delta_m(int J),
       Average_r_i(int J),Average_e_i(int J),
       Interpolate_U(double f[], double x1),Interpolate_V(double f[], double x1),
       Secret(double x),Deplet(double x1, double x2),GAUSS(),Deplet_V(int J),RAN(),
       VfrontPosition(double vf),Sample_from_distribution(double,int,double,double,double,double),
       LogNormal_distr(double,double,double,double),Normal_distr(double,double,double,double),
       Uniform_distr(double,double);
int CountCells(double x1, double x2),SIGN(double z),CellsInTheClot();
double dt,TIME,dX_u,dX_v,k_S,
       alpha,kappa,D_v,S_tax,alpha_0,kappa_0,D_v_0,S_tax_0,k_v,k_u;
int Ncells_n,Ncells_n_1,Ncells_max,Ncells_0,DISTR_rmax,n_runs,n_0,i_run=1;
short STOP=0;
char coeff[50], directory[50], inputfile[50];
double Umin=0.,Umax=10.,DMmin=0.,DMmax=0.1,Vmin=0.,Vmax=10.,Cmin=0.,Cmax=1.,
       PenDepth=0.,CellsInvaded=0.;

double Rmin=0.,Rmax=4.,rmax_0,MEAN_rmax,SIGMA_rmax,
       Pmin=0.,Pmax=4.,pmax_0,MEAN_pmax,SIGMA_pmax,
       Emin=0.,Emax=4.,emax_0,MEAN_emax,SIGMA_emax;
int DISTR_rmax,DISTR_pmax,DISTR_emax,same_rmax,same_pmax,same_emax,Model;

// ***************************************************************************************
// *** INITIALIZATION OF MODEL PARAMETRS, INITIAL CONDITIONS *****************************
void ModelParameters(){
int i,j;
char filename[50],variable[50];
FILE *init,*distr;

   InitialRandomSeed();

   if(i_run==1){
      init=fopen(inputfile,"r");
      fscanf(init,"%s %d",&variable, &n_0);
      fscanf(init,"%s %lf",&variable, &k_S);

      fscanf(init,"%s %lf",&variable, &alpha_0);
      fscanf(init,"%s %lf",&variable, &kappa_0);
      fscanf(init,"%s %lf",&variable, &D_v_0);	// mm^2/h
      fscanf(init,"%s %lf",&variable, &S_tax_0);// mm/h
      fscanf(init,"%s %lf",&variable, &k_v);	// 1/h
      fscanf(init,"%s %lf",&variable, &k_u);

      fscanf(init,"%s %d",&variable,  &Model);

      fscanf(init,"%s %lf",&variable, &MEAN_rmax);// used for both Normal and LogNormal distributions
      fscanf(init,"%s %lf",&variable, &SIGMA_rmax);// used for both Normal and LogNormal distributions
      fscanf(init,"%s %d" ,&variable, &DISTR_rmax);// 0: Constant; 1: Normal; 2: LogNormal; 3: Uniform
      fscanf(init,"%s %lf",&variable, &rmax_0);// if rmax is constant

      fscanf(init,"%s %lf",&variable, &MEAN_pmax);
      fscanf(init,"%s %lf",&variable, &SIGMA_pmax);
      fscanf(init,"%s %d" ,&variable, &DISTR_pmax);
      fscanf(init,"%s %lf",&variable, &pmax_0);

      fscanf(init,"%s %lf",&variable, &MEAN_emax);
      fscanf(init,"%s %lf",&variable, &SIGMA_emax);
      fscanf(init,"%s %d" ,&variable, &DISTR_emax);
      fscanf(init,"%s %lf",&variable, &emax_0);

      fscanf(init,"%s %d",&variable, &same_rmax);
      fscanf(init,"%s %d",&variable, &same_pmax);
      fscanf(init,"%s %d",&variable, &same_emax);

      fscanf(init,"%s %d",&variable, &n_runs);
      fclose(init);
   }

   dX_u=Xmax/((float)(N_U));
   dX_v=Xmax/((float)(N_V));
   Ncells_0=n_0*x_dermis*dX_u;	// the initial number of cells in the dermis region
   Ncells_n=Ncells_0;
   Ncells_n_1=Ncells_n;
   Ncells_max=100*Ncells_0;
   dt=0.1;

   Create_Arrays();

   TIME=0.;
   for(j=0;j<N_U;j++){
      X_u[j]=(0.5+j)*dX_u;
      U_n[j]=0.;
      U_n_1[j]=U_n[j];
      gradU_n[j]=0.;
      gradU_n_1[j]=gradU_n[j];
   }
   AU[0]=0.;
   BU[0]=1.;
   CU[0]=0.;
   DU[0]=0.;
   for(j=1;j<N_U-1;j++){
      AU[j]=-D_u/(dX_u*dX_u);
      BU[j]=(1./dt)+2.*D_u/(dX_u*dX_u);
      CU[j]=AU[j];
   }
   AU[N_U-1]=-2.*D_u/(dX_u*dX_u);
   BU[N_U-1]=(1./dt)+2.*D_u/(dX_u*dX_u);
   CU[N_U-1]=0.;

   for(j=0;j<Ncells_n;j++){
      (CELL_n+j)->x = x_dermis*RAN();
      (CELL_n+j)->y = RAN();
      // Initial distribution is defined in the input file
      (CELL_n+j)->rmax = 
         Sample_from_distribution(rmax_0,DISTR_rmax,MEAN_rmax,SIGMA_rmax,Rmin,Rmax);
      (CELL_n+j)->pmax = 
         Sample_from_distribution(pmax_0,DISTR_pmax,MEAN_pmax,SIGMA_pmax,Pmin,Pmax);
      (CELL_n+j)->emax = 
         Sample_from_distribution(emax_0,DISTR_emax,MEAN_emax,SIGMA_emax,Emin,Emax);

      (CELL_n+j)->r_i = 0.;
      (CELL_n+j)->e_i = 0.;
      (CELL_n+j)->delta_m = 0.;
      (CELL_n+j)->nabla_m = 0.;
      (CELL_n_1+j)->x = (CELL_n+j)->x;
      (CELL_n_1+j)->y = (CELL_n+j)->y;
      (CELL_n_1+j)->rmax = (CELL_n+j)->rmax;
      (CELL_n_1+j)->pmax = (CELL_n+j)->pmax;
      (CELL_n_1+j)->emax = (CELL_n+j)->emax;
      (CELL_n_1+j)->r_i = (CELL_n+j)->r_i;
      (CELL_n_1+j)->e_i = (CELL_n+j)->e_i;
      (CELL_n_1+j)->delta_m = (CELL_n+j)->delta_m;
      (CELL_n_1+j)->nabla_m = (CELL_n+j)->nabla_m;
   }
   for(j=0;j<N_V;j++)
      X_v[j]=(0.5+j)*dX_v;

   sprintf(filename,"%s/rmax.distr",coeff);
   distr=fopen(filename,"w");
   for(j=0;j<Ncells_n;j++)
      fprintf(distr,"\n %f",(CELL_n+j)->rmax);
   fclose(distr);

   sprintf(filename,"%s/pmax.distr",coeff);
   distr=fopen(filename,"w");
   for(j=0;j<Ncells_n;j++)
      fprintf(distr,"\n %f",(CELL_n+j)->pmax);
   fclose(distr);

   sprintf(filename,"%s/emax.distr",coeff);
   distr=fopen(filename,"w");
   for(j=0;j<Ncells_n;j++)
      fprintf(distr,"\n %f",(CELL_n+j)->emax);
   fclose(distr);

   CellDensity();
   CellDensity();
}

double Sample_from_distribution(double r0, int d, double m, double s, double x1, double x2){
double r;
   switch (d){
      case 0:
         r=r0;
         break;
      case 1:
         r=Normal_distr(m,s,x1,x2);
         break;
      case 2:
         r=LogNormal_distr(m,s,x1,x2);
         break;
      case 3:
         r=Uniform_distr(x1,x2);
         break;
   }
   return r;
}

void Create_Arrays(){

   X_u	    =(ARRAY) calloc(N_U,sizeof(TYPE));
   X_v      =(ARRAY) calloc(N_U,sizeof(TYPE));
   U_n	    =(ARRAY) calloc(N_U,sizeof(TYPE));
   U_n_1    =(ARRAY) calloc(N_U,sizeof(TYPE));
   gradU_n  =(ARRAY) calloc(N_U,sizeof(TYPE));
   gradU_n_1=(ARRAY) calloc(N_U,sizeof(TYPE));
   V_n	    =(ARRAY) calloc(N_V,sizeof(TYPE));
   V_n_1    =(ARRAY) calloc(N_V,sizeof(TYPE));
   AU	    =(ARRAY) calloc(N_U,sizeof(TYPE));
   BU	    =(ARRAY) calloc(N_U,sizeof(TYPE));
   CU	    =(ARRAY) calloc(N_U,sizeof(TYPE));
   DU	    =(ARRAY) calloc(N_U,sizeof(TYPE));

   KILL     =(INDEX) calloc(Ncells_max,sizeof(TYPE1));
   DIVIDE   =(INDEX) calloc(Ncells_max,sizeof(TYPE1));

   CELL_n   =(struct cells *)malloc((unsigned)(Ncells_max*9*sizeof(double)));
   CELL_n_1 =(struct cells *)malloc((unsigned)(Ncells_max*9*sizeof(double)));
   temp     =(struct TEMP *) malloc((unsigned)(Ncells_max*5*sizeof(double)));
}


// ***************************************************************************************
// *** STORE DATA AT THE PREVIOUS STEP n-1 AND SOLVE LIGAND PDE FOR U VALUES AT STEP n ***
void SolvePDE(){
int j;

   for(j=0;j<N_U;j++)
      U_n_1[j]=U_n[j];
   for(j=1;j<N_U;j++){
      DU[j]= k_S*Secret(X_u[j]) - k_u*U_n_1[j] - 
      		k_v*Deplet(X_u[j]-0.5*dX_u,X_u[j]+0.5*dX_u)/(n_0*dX_u*dX_u) + (U_n_1[j]/dt);
   }
   Tridiag(AU,BU,CU,U_n,DU);
   ComputeGradient(U_n_1,gradU_n_1);
   ComputeGradient(U_n,gradU_n);
}
// ***************************************************************************************
// *** STORE DATA AT THE PREVIOUS STEP n-1 AND UPDATE CELLS AT STEP n ********************
void MoveCells(){
double C,f,u_i,gu_i,a,P1,P2,P,k_mitosis,k_death,v,m_i,y,e_i,r_i,e_max,r_max,p_max,q0,q1,q2;
int i,N,M,k,j,n=0,m=0;

   Ncells_n_1=Ncells_n;
   for(i=0;i<Ncells_n_1;i++){
      (CELL_n_1+i)->x = (CELL_n+i)->x;
      (CELL_n_1+i)->y = (CELL_n+i)->y;
      (CELL_n_1+i)->rmax = (CELL_n+i)->rmax;
      (CELL_n_1+i)->pmax = (CELL_n+i)->pmax;
      (CELL_n_1+i)->emax = (CELL_n+i)->emax;
      (CELL_n_1+i)->e_i = (CELL_n+i)->e_i;
      (CELL_n_1+i)->r_i = (CELL_n+i)->r_i;
      (CELL_n_1+i)->delta_m = (CELL_n+i)->delta_m;
      (CELL_n_1+i)->nabla_m = (CELL_n+i)->nabla_m;
   }

   for(i=0;i<Ncells_n_1;i++){
      KILL[i]=-1;
      DIVIDE[i]=-1;
   }

   for(i=0;i<Ncells_n_1;i++){
      v = Interpolate_V(V_n_1,(CELL_n_1+i)->x);
      if( ((CELL_n_1+i)->x) < x_dermis ){
         k_death=k_d*(1.-V_0/v);
	 if(k_death<0.)
	    k_death=0.;
      }else
         k_death=k_d;
      P1=1.-exp(-k_death*dt);
      k_mitosis=mu_m*((CELL_n_1+i)->r_i)*
         (1.-(v/V_star)*(v/V_star)*(v/V_star))/(GAMMA+((CELL_n_1+i)->r_i));
      P2=1.-exp(-k_mitosis*dt);
      P=P1+P2;

      if(P>=1.)
         printf("\n P=%f",P);

      if(P>RAN()){
         if( (P1/P)>RAN() ){
	    KILL[n]=i;
	    n++;
	 }else{
	    DIVIDE[m]=i;
	    m++;
	 }
      }
      else{
         D_v=D_v_0;
         S_tax=S_tax_0;
         m_i=((CELL_n_1+i)->e_i + m0*((CELL_n_1+i)->emax - (CELL_n_1+i)->e_i))/(CELL_n_1+i)->pmax;
	 C=SIGN(RAN()-0.5)*sqrt(2.*D_v*m_i/dt)+
            (S_tax*((CELL_n_1+i)->delta_m)-D_v*((CELL_n_1+i)->nabla_m));
	 (CELL_n+i)->x = (CELL_n_1+i)->x + C*dt;
         if( ((CELL_n+i)->x)<0. )
            (CELL_n+i)->x = -(CELL_n+i)->x;
         if( ((CELL_n+i)->x)>Xmax )
            (CELL_n+i)->x = Xmax - ((CELL_n+i)->x - Xmax);
      }
   }

   if(n>0 || m>0){
      Ncells_n=Ncells_n_1-n+m;
      N=0; M=0; k=0;
      for(j=0;j<Ncells_n_1;j++){
   	 if(j!=KILL[N]){
   	    (temp+k)->x = (CELL_n+j)->x;
   	    (temp+k)->y = (CELL_n+j)->y;
   	    (temp+k)->rmax = (CELL_n+j)->rmax;	// parent cell
   	    (temp+k)->pmax = (CELL_n+j)->pmax;	// parent cell
   	    (temp+k)->emax = (CELL_n+j)->emax;	// parent cell
     	    if(j==DIVIDE[M]){
   	      (temp+(Ncells_n_1-n)+M)->x = (CELL_n+j)->x;
   	      (temp+(Ncells_n_1-n)+M)->y = (CELL_n+j)->y;

              if(same_rmax==1){	// *** child cells with the same rmax
                 f=(temp+k)->rmax;
   	             (temp+(Ncells_n_1-n)+M)->rmax = f;
   	             (temp+k)->rmax = f;
              }else{
	             q0=(temp+k)->rmax;	// rmax of the parent cell
				 switch (Model){
                    case 1:	// **** Model I *****************************
				       q1=Sample_from_distribution(rmax_0,DISTR_rmax,1.,SIGMA_rmax,Rmin,Rmax);
                       q2=Sample_from_distribution(rmax_0,DISTR_rmax,1.,SIGMA_rmax,Rmin,Rmax);
					break;
					case 2:	// **** Model II *****************************
				       q1=Sample_from_distribution(rmax_0,DISTR_rmax,q0,SIGMA_rmax,Rmin,Rmax);
                       q2=Sample_from_distribution(rmax_0,DISTR_rmax,q0,SIGMA_rmax,Rmin,Rmax);
					break;
                    case 3:	// **** Model III *****************************
                       q1 = q0;
   	                   q2 = q0;
					break;
			     }
				 (temp+k)->rmax = q1;
				 (temp+(Ncells_n_1-n)+M)->rmax = q2;
              }

              if(same_pmax==1){	// *** child cells with the same pmax
                 f=(temp+k)->pmax;
   	             (temp+(Ncells_n_1-n)+M)->pmax = f;
   	             (temp+k)->pmax = f;
              }else{
	             q0=(temp+k)->pmax;	// pmax of the parent cell
				 switch (Model){
                    case 1:	// **** Model I *****************************
                       q1=Sample_from_distribution(pmax_0,DISTR_pmax,1.,SIGMA_pmax,Pmin,Pmax);
                       q2=Sample_from_distribution(pmax_0,DISTR_pmax,1.,SIGMA_pmax,Pmin,Pmax);
					break;
					case 2:	// **** Model II *****************************
                       q1=Sample_from_distribution(pmax_0,DISTR_pmax,q0,SIGMA_pmax,Pmin,Pmax);
                       q2=Sample_from_distribution(pmax_0,DISTR_pmax,q0,SIGMA_pmax,Pmin,Pmax);
					break;
					case 3:	// **** Model III *****************************
					   q1 = q0;
					   q2 = q0;
					break;
                 }
                 (temp+k)->pmax = q1;
   	             (temp+(Ncells_n_1-n)+M)->pmax = q2;
              }

              if(same_emax==1){	// *** child cells with the same emax
                 f=(temp+k)->emax;
   	             (temp+(Ncells_n_1-n)+M)->emax = f;
   	             (temp+k)->emax = f;
              }else{		// *** change emax to random for both new cells
	             q0=(temp+k)->emax;	// emax of the parent cell
				 switch (Model){
				    case 1:	// **** Model I *****************************
                       q1=Sample_from_distribution(emax_0,DISTR_emax,1.,SIGMA_emax,Emin,Emax);
                       q2=Sample_from_distribution(emax_0,DISTR_emax,1.,SIGMA_emax,Emin,Emax);
					break;
                    case 2:	// **** Model II *****************************
                       q1=Sample_from_distribution(emax_0,DISTR_emax,q0,SIGMA_emax,Emin,Emax);
                       q2=Sample_from_distribution(emax_0,DISTR_emax,q0,SIGMA_emax,Emin,Emax);
					break;
				    case 3:	// **** Model III *****************************
					   q1 = q0;
					   q2 = q0;
					break;
                 }
                 (temp+k)->emax = q1;
   	             (temp+(Ncells_n_1-n)+M)->emax = q2;
              }

   	      y = (temp+(Ncells_n_1-n)+M)->y + 0.01;
	      if( y<1. )
	         (temp+(Ncells_n_1-n)+M)->y = y;
     	      y = ((temp+k)->y) - 0.01;
              if( y>0. )
	         (temp+k)->y = y;

	      if(M<m-1)
     		 M++;
     	    }
	    k++;
   	 }else{
     	    if(N<n-1)
     	       N++;
   	 }
      }

      for(i=0;i<Ncells_n;i++){
         (CELL_n+i)->x = (temp+i)->x;
         (CELL_n+i)->y = (temp+i)->y;
         (CELL_n+i)->rmax = (temp+i)->rmax;
         (CELL_n+i)->pmax = (temp+i)->pmax;
         (CELL_n+i)->emax = (temp+i)->emax;
      }
   }

   for(i=0;i<Ncells_n;i++){
      alpha=alpha_0;
      kappa=kappa_0;
      u_i = Interpolate_U(U_n, (CELL_n+i)->x);
      gu_i = Interpolate_U(gradU_n, (CELL_n+i)->x);
      f = (2.*u_i + u_i*u_i)*gu_i/((1.+u_i + u_i*u_i)*(1.+u_i + u_i*u_i));
      (CELL_n+i)->r_i = ReceptorLevel(U_n, (CELL_n+i)->x, (CELL_n+i)->rmax);

      r_i=(CELL_n+i)->r_i;
      e_max=(CELL_n+i)->emax;
      a=e_max+kappa+alpha*r_i;
      e_i=0.5*(a - sqrt(-4.*alpha*r_i*e_max + a*a));
      r_max=(CELL_n+i)->rmax;
      p_max=(CELL_n+i)->pmax;

      (CELL_n+i)->e_i = e_i;
      (CELL_n+i)->nabla_m = ((1.-m0)/p_max)*alpha*r_max*(e_max-e_i)*f/(a-2.*e_i);
      (CELL_n+i)->delta_m = (0.5*alpha*r_max/p_max)*(e_max-e_i)*epsilon*f/(kappa+e_max-e_i);
   }
}

// ***************************************************************************************
// *** STORE DATA AT THE PREVIOUS STEP n-1 AND UPDATE CELL DENSITY AT STEP n *************
// *** V IS COMPUTED ON A ROUGH GRID TO SMOOTH OUT FLUCTUATIONS **************************
void CellDensity(){
int i;

   for(i=0;i<N_V;i++)
      V_n_1[i]=V_n[i];

   for(i=0;i<N_V;i++)
      V_n[i]=CountCells(X_v[i]-0.5*dX_v,X_v[i]+0.5*dX_v)/(dX_u*dX_v)/n_0;
}
int CountCells(double x1, double x2){
int i,j=0;

   for(i=0;i<Ncells_n;i++){
      if( ((CELL_n+i)->x)>=x1 && ((CELL_n+i)->x)<x2 )
         j++;
   }
   return j;
}

// ***************************************************************************************
// *** THIS FUNCTION IS USED AT n-1 STEP ONLY ********************************************
double Deplet(double x1, double x2){
double R=0.;
int I;

   for(I=0;I<Ncells_n_1;I++){
      if( ((CELL_n_1+I)->x)>=x1 && ((CELL_n_1+I)->x)<x2 )
         R=R+ReceptorLevel(U_n_1,(CELL_n_1+I)->x,(CELL_n_1+I)->rmax);
   }
   return R;
}

// ***************************************************************************************
// *** THESE FUNCTIONS ARE USED AT BOTH n AND n-1 STEPS **********************************
double ReceptorLevel(double U[], double x_i, double rmax_i){
double r, u_i;
   u_i = Interpolate_U(U, x_i);
   r = rmax_i*u_i*u_i / (1. + u_i + u_i*u_i);
   return r;
}
void ComputeGradient(double U[],double gU[]){
int j;

   gU[0] = (U[1]-U[0])/dX_u;
   for(j=1;j<N_U-1;j++)
      gU[j] = (U[j+1]-U[j-1])/(2.*dX_u);
   gU[N_U-1] = (U[N_U-1]-U[N_U-2])/dX_u;
}
double Interpolate_U(double f[], double x){
double F;
int j;

   j=(int)(x/dX_u - 0.5);   
   F=((f[j+1]-f[j])/(X_u[j+1]-X_u[j]))*(x-X_u[j])+f[j];
   return F;
}
double Interpolate_V(double f[], double x){
double F;
int j;

   j=(int)(x/dX_v - 0.5);   
   F=((f[j+1]-f[j])/(X_v[j+1]-X_v[j]))*(x-X_v[j])+f[j];
   return F;
}
double Secret(double x){
   if(x<x_dermis)
      return 0.;
   else
      return 1.;
}

int SIGN(double z){
   if(z<0.)
      return -1;
   else
      return 1;
}
double Tridiag(double a[],double b[],double c[],double f[],double d[])
{
   int k;
   double gam[N_U],bet;
   if(b[0]==0){
     printf("\n  Rewrite Equations !");
     exit(0);
   }
   bet=b[0];
   f[0]=d[0]/bet;
   for(k=1;k<N_U;k++){
     gam[k]=c[k-1]/bet;
     bet=b[k]-a[k]*gam[k];
     if(bet==0){
       printf("\n  Thomas Failed !");
       exit(0);
     }
     f[k]=(d[k]-a[k]*f[k-1])/bet;
   }
   for(k=N_U-2;k>=0;k--)
     f[k]=f[k]-gam[k+1]*f[k+1];
}
double GAUSS(){
static int iset = 0;
static double gset;
double  fac, v1, v2, rsquare=2.0;

   if(iset == 0){
      while(rsquare>1.0 || rsquare==0.0){
        v1 = 2*RAN()-1.0;
        v2 = 2*RAN()-1.0;
        rsquare = v1*v1 + v2*v2;
      }; 
      fac=sqrt(-2.0*log(rsquare)/rsquare);
      gset=v1*fac;
      iset=1;
      return v2*fac;
   }else{
      iset=0;
      return gset;
   }
}
double RAN(){
   return drand48();
}

double Normal_distr(double m, double s, double x1, double x2){
double P;
   do{
      //P = m+s*GAUSS();
      P = m*(1.+s*GAUSS());
   }while( P<x1 || P>x2 );
   return P;
}
double LogNormal_distr(double m, double s, double x1, double x2){
double P;
   do{
      P = m*exp(s*GAUSS());
   }while( P<x1 || P>x2 );
   return P;
}
double Uniform_distr(double x1, double x2){
   return x1+(x2-x1)*RAN();
}

double Average_r_max(int J){
double x, r=0., x1=(X_v[J]-0.5*dX_v), x2=(X_v[J]+0.5*dX_v);
int k,n=0;

   for(k=0;k<Ncells_n;k++){
      x=(CELL_n+k)->x;
      if( x>=x1 && x<x2 ){
         r=r+(CELL_n+k)->rmax;
	 n++;
      }
   }
   if( n>0 )
      return r/n;
   else
      return 0.;
}
double Average_p_max(int J){
double x, r=0., x1=(X_v[J]-0.5*dX_v), x2=(X_v[J]+0.5*dX_v);
int k,n=0;

   for(k=0;k<Ncells_n;k++){
      x=(CELL_n+k)->x;
      if( x>=x1 && x<x2 ){
         r=r+(CELL_n+k)->pmax;
	 n++;
      }
   }
   if( n>0 )
      return r/n;
   else
      return 0.;
}
double Average_e_max(int J){
double x, r=0., x1=(X_v[J]-0.5*dX_v), x2=(X_v[J]+0.5*dX_v);
int k,n=0;

   for(k=0;k<Ncells_n;k++){
      x=(CELL_n+k)->x;
      if( x>=x1 && x<x2 ){
         r=r+(CELL_n+k)->emax;
	 n++;
      }
   }
   if( n>0 )
      return r/n;
   else
      return 0.;
}
double Average_Delta_m(int J){
double x, r=0., x1=(X_v[J]-0.5*dX_v), x2=(X_v[J]+0.5*dX_v);
int k,n=0;

   for(k=0;k<Ncells_n;k++){
      x=(CELL_n+k)->x;
      if( x>=x1 && x<x2 ){
         r=r+(CELL_n+k)->delta_m;
	 n++;
      }
   }
   if( n>0 )
      return r/n;
   else
      return 0.;
}
double Average_r_i(int J){
double x, r=0., x1=(X_v[J]-0.5*dX_v), x2=(X_v[J]+0.5*dX_v);
int k,n=0;

   for(k=0;k<Ncells_n;k++){
      x=(CELL_n+k)->x;
      if( x>=x1 && x<x2 ){
         r=r+(CELL_n+k)->r_i;
	 n++;
      }
   }
   if( n>0 )
      return r/n;
   else
      return 0.;
}
double Average_e_i(int J){
double x, r=0., x1=(X_v[J]-0.5*dX_v), x2=(X_v[J]+0.5*dX_v);
int k,n=0;

   for(k=0;k<Ncells_n;k++){
      x=(CELL_n+k)->x;
      if( x>=x1 && x<x2 ){
         r=r+(CELL_n+k)->e_i;
	 n++;
      }
   }
   if( n>0 )
      return r/n;
   else
      return 0.;
}


void InitialRandomSeed(){
      RAN();
}

// ***************************************************************************************
// *** DEFINES x-POSITION OF THE v-FRONT AT v=vf *****************************************
double VfrontPosition(double vf){
int i;
double x;
   for(i=N_V-1;i>=1;i--){
      if( V_n[i]<vf && V_n[i-1]>vf ){
	 x=(vf-V_n[i-1])*(X_v[i]-X_v[i-1])/(V_n[i]-V_n[i-1])+X_v[i-1];
	 break;
      }
   }
   return x;
}

int CellsInTheClot(){
int i,n=0;

   for(i=0;i<Ncells_n;i++){
      if( ((CELL_n+i)->x)>x_dermis )
         n++;
   }
   return n;
}


// **********************************************************************************************
// **********************************************************************************************
// **********************************************************************************************
// **********************************************************************************************
// *********** GUI parameters and variables *****************************************************
// **********************************************************************************************
// **********************************************************************************************
// **********************************************************************************************
// **********************************************************************************************
#define DIMENSION_x	650
#define DIMENSION_y	220
Boolean Loop(XtPointer client_data);
Widget toplevel,iter_box;
GC context;
Window window_U,window_V,window_C;
Display *display;
XColor color[256];
static XtAppContext app_con;

XPoint *PDGF,*DENS;
int  XtoGaphics(double X),UtoGaphics(double U),DMtoGaphics(double U),VtoGaphics(double V),
     CtoGaphics(double C),RtoGaphics(double R);
void Quit(Widget w, XtPointer call_data, XtPointer client_data),
     Redraw(Widget w, XtPointer call_data, XtPointer client_data),
     Start(Widget w,XtPointer client_data,XtPointer call_data),
     Calc(Widget w,XtPointer client_data, XtPointer call_data),
     Show(Widget w,XtPointer client_data, XtPointer call_data),
     Print_Data(Widget w, XtPointer call_data, XtPointer client_data),
     DrawFrame(int,int),ColorScale(),GraphicalParameters(),PlotData(int,int),Print_Time(),
     CheckMaxLevels();
Boolean Loop(XtPointer client_data);
short on_off=1,show=0;
int GDx,GDy,GLx,GLy,print_out;

int main(int argc, char **argv){
Widget main_box,central_box,plots_box_U,plots_box_V,plots_box_C,start,quit,image,calc,redraw,print;
Arg args[10];
int i;
static String display_on_off[]={"Image Off"," Image On",NULL};
static String calc_on_off[]={"Calc Off","  Calc On",NULL};

   sprintf(inputfile,"%s",argv[1]);
   sprintf(directory,"%s",argv[2]);
   mkdir(directory,0777);
   sprintf(coeff,"%s",directory);
   mkdir(coeff,0777);
   ModelParameters();

   toplevel = XtAppInitialize(&app_con,"toplevel",NULL,ZERO,&argc,argv,NULL,NULL,ZERO);
   i=0;
   XtSetArg(args[i],XtNx,0); i++;
   XtSetArg(args[i],XtNy,0); i++;
   XtSetValues(toplevel,args,i);
   i=0;
   XtSetArg(args[i],XtNx,0); i++;
   XtSetArg(args[i],XtNy,0); i++;
   main_box=XtCreateManagedWidget("main_box", boxWidgetClass,toplevel,args,i);
   central_box=XtCreateManagedWidget("central_box", boxWidgetClass,main_box,args,0);
   i=0;
   XtSetArg(args[i],XtNlist,display_on_off); i++;
   XtSetArg(args[i],XtNwidth,200); i++;
   image=XtCreateManagedWidget("display_on_off",listWidgetClass,central_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlist,calc_on_off); i++;
   XtSetArg(args[i],XtNwidth,200); i++;
   calc=XtCreateManagedWidget("calc_on_off",listWidgetClass,central_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlabel,"Redraw"); i++;
   redraw=XtCreateManagedWidget("redraw",commandWidgetClass,central_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlabel,"Start"); i++;
   start=XtCreateManagedWidget("start",commandWidgetClass,central_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlabel,"Quit"); i++;
   quit=XtCreateManagedWidget("quit",commandWidgetClass,central_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlabel,"time [s]: 0.0"); i++;
   XtSetArg(args[i],XtNwidth,120); i++;
   iter_box=XtCreateManagedWidget(NULL,labelWidgetClass,central_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlabel,"Position"); i++;
   print=XtCreateManagedWidget("print",commandWidgetClass,central_box,args,i);

   i=0;
   XtSetArg(args[i],XtNwidth,DIMENSION_x); i++;
   XtSetArg(args[i],XtNheight,DIMENSION_y); i++;
   plots_box_U=XtCreateManagedWidget("plots_box_U",coreWidgetClass,central_box,args,i);
   i=0;
   XtSetArg(args[i],XtNwidth,DIMENSION_x); i++;
   XtSetArg(args[i],XtNheight,DIMENSION_y); i++;
   plots_box_V=XtCreateManagedWidget("plots_box_V",coreWidgetClass,central_box,args,i);
   i=0;
   XtSetArg(args[i],XtNwidth,DIMENSION_x); i++;
   XtSetArg(args[i],XtNheight,DIMENSION_y); i++;
   plots_box_C=XtCreateManagedWidget("plots_box_C",coreWidgetClass,central_box,args,i);

   XtAddCallback(start,XtNcallback,Start,NULL);
   XtAddCallback(quit,XtNcallback,Quit,NULL);
   XtAddCallback(redraw,XtNcallback,Redraw,NULL);
   XtAddCallback(calc,XtNcallback,Calc,NULL);
   XtAddCallback(image,XtNcallback,Show,NULL);
   XtAddCallback(print,XtNcallback,Print_Data,NULL);

   XtRealizeWidget(toplevel);
   
   display=XtDisplay(plots_box_U);
   window_U=XtWindow(plots_box_U);
   context=XCreateGC(display,window_U,0,NULL);

   window_V=XtWindow(plots_box_V);
   window_C=XtWindow(plots_box_C);

   ColorScale();
   GraphicalParameters();
   Start(NULL,NULL,NULL);

   XtAppMainLoop(app_con);
}


void Start(Widget w, XtPointer client_data, XtPointer call_data){
   ModelParameters();
   print_out=0;
   Redraw(NULL,NULL,NULL);
   XtAppAddWorkProc(app_con,Loop,NULL);
}
void Calc(Widget w,XtPointer client_data, XtPointer call_data){
   XawListReturnStruct *item = (XawListReturnStruct*)call_data;
   on_off=item->list_index;
}
void Show(Widget w,XtPointer client_data, XtPointer call_data){
   XawListReturnStruct *item = (XawListReturnStruct*)call_data;
   show=item->list_index;
   Redraw(NULL,NULL,NULL);
}
void Quit(Widget w, XtPointer call_data, XtPointer client_data){
   Average_and_print_depth();
   exit(0);
}
void Redraw(Widget w, XtPointer call_data, XtPointer client_data){
   CheckMaxLevels();
   DrawFrame(Rmin,Rmax);
   PlotData(Rmin,Rmax);
}
Boolean Loop(XtPointer client_data){
int n;
   switch (on_off){
      case 0:
      break;
      case 1:
	SolvePDE();
	MoveCells();
	CellDensity();
	TIME=TIME+dt;
        print_out++;
	if( print_out==10 ){
	    Print_Time();
	    print_out=0;
	    //PrintMovieData();
	}
	if(show==1){
	   CheckMaxLevels();
	   PlotData(Rmin,Rmax);
	}
	if(TIME>=TIME_final){
	   Print_Data(NULL,NULL,NULL);
           if(i_run==n_runs)
	      Quit(NULL,NULL,NULL);
           i_run++;
           ModelParameters();
	}
      break;
   }
   return(False);
}

void GraphicalParameters(){
int j;
   PDGF=(XPoint *) malloc((unsigned) 2*(N_U)*sizeof(short));
   DENS=(XPoint *) malloc((unsigned) 2*(N_U)*sizeof(short));
   GDx=(int)(0.07*DIMENSION_x);
   GDy=(int)(0.11*DIMENSION_y);
   GLx=(int)(0.85*DIMENSION_x);
   GLy=(int)(0.8*DIMENSION_y);
   for(j=0;j<N_U;j++)
      (PDGF+j)->x = XtoGaphics(X_u[j]);
   for(j=0;j<N_V;j++)
      (DENS+j)->x = XtoGaphics(X_v[j]);
}
void DrawFrame(int ymin, int ymax){
char Label[50];
int i,chub=4;

// *** DRAW FRAME FOR LIGAND COMCENTRATION, U ************************************************************************
   XSetForeground(display,context,color[196].pixel);
   XFillRectangle(display,window_U,context,0,0,DIMENSION_x,DIMENSION_y);
   XSetForeground(display,context,1);
   XDrawRectangle(display,window_U,context,GDx-1,GDy-1,GLx+2,GLy+2);
   for(i=0;i<=10;i++){
      XDrawLine(display,window_U,context,
      	XtoGaphics(0.1*i*Xmax),GDy+GLy+chub,XtoGaphics(0.1*i*Xmax),GDy+GLy+1);
      XDrawLine(display,window_U,context,
      	GDx-chub,UtoGaphics(0.1*i*(Umax-Umin)+Umin),GDx-1,UtoGaphics(0.1*i*(Umax-Umin)+Umin));
      XDrawLine(display,window_U,context,
      	GDx+GLx+chub,UtoGaphics(0.1*i*(Umax-Umin)+Umin),GDx+GLx+1,UtoGaphics(0.1*i*(Umax-Umin)+Umin));
   }
   XSetForeground(display,context,1);
   for(i=0;i<=10;i++){
      sprintf(Label,"%.2f",(0.1*i*Xmax));
      XDrawString(display,window_U,context,XtoGaphics(0.1*i*Xmax)-10,GDy+GLy+16,Label,strlen(Label));
   }
   XSetForeground(display,context,color[180].pixel);
   for(i=0;i<=10;i++){
      sprintf(Label,"%.2f",(0.1*i*(Umax-Umin)+Umin));
      XDrawString(display,window_U,context,GDx+GLx+10,UtoGaphics(0.1*i*(Umax-Umin)+Umin)+5,Label,strlen(Label));
   }
   sprintf(Label,"PDGF concentration, u");
   XDrawString(display,window_U,context,GDx+GLx-130,GDy-8,Label,strlen(Label));
   XSetForeground(display,context,1);
   for(i=0;i<=10;i++){
      sprintf(Label,"%.2f",(0.1*i*(DMmax-DMmin)+DMmin));
      XDrawString(display,window_U,context,GDx-35,DMtoGaphics(0.1*i*(DMmax-DMmin)+DMmin)+5,Label,strlen(Label));
   }
   sprintf(Label,"Chemotactic signaling at cells, delta m");
   XDrawString(display,window_U,context,GDx,GDy-8,Label,strlen(Label));
// *** DRAW FRAME FOR CELL DENSITY, V ********************************************************************************
   XSetForeground(display,context,color[196].pixel);
   XFillRectangle(display,window_V,context,0,0,DIMENSION_x,DIMENSION_y);
   XSetForeground(display,context,1);
   XDrawRectangle(display,window_V,context,GDx-1,GDy-1,GLx+2,GLy+2);
   for(i=0;i<=10;i++){
      XDrawLine(display,window_V,context,
      	XtoGaphics(0.1*i*Xmax),GDy+GLy+chub,XtoGaphics(0.1*i*Xmax),GDy+GLy+1);
      XDrawLine(display,window_V,context,
      	GDx-chub,VtoGaphics(0.1*i*(Vmax-Vmin)+Vmin),GDx-1,VtoGaphics(0.1*i*(Vmax-Vmin)+Vmin));
      XDrawLine(display,window_V,context,
      	GDx+GLx+chub,RtoGaphics(0.1*i*(ymax-ymin)+ymin),GDx+GLx+1,
        RtoGaphics(0.1*i*(ymax-ymin)+ymin));
   }
   XSetForeground(display,context,1);
   for(i=0;i<=10;i++){
      sprintf(Label,"%.2f",(0.1*i*(Vmax-Vmin)+Vmin));
      XDrawString(display,window_V,context,GDx-35,VtoGaphics(0.1*i*(Vmax-Vmin)+Vmin)+5,Label,strlen(Label));
      sprintf(Label,"%.2f",0.1*i*Xmax);
      XDrawString(display,window_V,context,XtoGaphics(0.1*i*Xmax)-10,GDy+GLy+16,Label,strlen(Label));
   }
   sprintf(Label,"Fibroblast density, v");
   XDrawString(display,window_V,context,GDx,GDy-8,Label,strlen(Label));

   XSetForeground(display,context,color[180].pixel);
   sprintf(Label,"r max, averaged over cells");
   XDrawString(display,window_V,context,GDx+GLx-130,GDy-8,Label,strlen(Label));
   for(i=0;i<=10;i++){
      sprintf(Label,"%.2f",(0.1*i*(ymax-ymin)+ymin));
      XDrawString(display,window_V,context,GDx+GLx+10,
         RtoGaphics(0.1*i*(ymax-ymin)+ymin)+5,Label,strlen(Label));
   }
   XSetForeground(display,context,1);

// *** DRAW FRAME FOR CELL POSITIONS *********************************************************************************
   XSetForeground(display,context,color[196].pixel);
   XFillRectangle(display,window_C,context,0,0,DIMENSION_x,DIMENSION_y);
   XSetForeground(display,context,1);
   XDrawRectangle(display,window_C,context,GDx-1,GDy-1,GLx+2,GLy+2);
   for(i=0;i<=10;i++){
      XDrawLine(display,window_C,context,
      	XtoGaphics(0.1*i*Xmax),GDy+GLy+chub,XtoGaphics(0.1*i*Xmax),GDy+GLy+1);
   }
   XSetForeground(display,context,1);
   for(i=0;i<=10;i++){
      sprintf(Label,"%.2f",(0.1*i*Xmax));
      XDrawString(display,window_C,context,XtoGaphics(0.1*i*Xmax)-10,GDy+GLy+16,Label,strlen(Label));
   }
   sprintf(Label,"Cell positions (y-coordinate is arbitrary)");
   XDrawString(display,window_C,context,GDx,GDy-8,Label,strlen(Label));

   sprintf(Label,"Colorscale for maximum receptor levels:");
   XDrawString(display,window_C,context,350,GDy-8,Label,strlen(Label));
   XSetForeground(display,context,1);
   XFillRectangle(display,window_C,context,GDx+GLx+18,0,DIMENSION_x-(GDx+GLx+18),DIMENSION_y);
   for(i=0;i<=20;i++){
      sprintf(Label,"%.1f",(0.05*i*(ymax-ymin)+ymin));
      XSetForeground(display,context,color[(int)(195*(1.-0.05*i))].pixel);
      XDrawString(display,window_C,context,
         GDx+GLx+22,CtoGaphics(1.1*(0.05*i*(Cmax-Cmin)+Cmin))+10,Label,strlen(Label));
   }
}
int XtoGaphics(double X){
   return (int)(GDx+X/Xmax*GLx);
}
int UtoGaphics(double U){
   return (int)(GDy+(1.-(U-Umin)/(Umax-Umin))*GLy);
}
int DMtoGaphics(double D){
   return (int)(GDy+(1.-(D-DMmin)/(DMmax-DMmin))*GLy);
}
int VtoGaphics(double V){
   return (int)(GDy+(1.-(V-Vmin)/(Vmax-Vmin))*GLy);
}
int CtoGaphics(double C){
   return (int)(GDy+(1.-(C-Cmin)/(Cmax-Cmin))*GLy);
}
int RtoGaphics(double R){
   return (int)(GDy+(1.-(R-Rmin)/(Rmax-Rmin))*GLy);
}
void PlotData(int ymin, int ymax){
int i,j;
char Label[50];
XPoint *FIBR;

// *** PDGF CONCENTRATION u AND CHEMOTACTIC SIGNALING delta m ***********************************************
   for(j=0;j<N_U;j++)
      (PDGF+j)->y = UtoGaphics(U_n_1[j]);
   XSetForeground(display,context,color[196].pixel);
   XDrawLines(display,window_U,context,PDGF,j,0);
   for(j=0;j<N_U;j++)
      (PDGF+j)->y = UtoGaphics(U_n[j]);
   XSetForeground(display,context,color[180].pixel);
   XDrawLines(display,window_U,context,PDGF,j,0);
   FIBR=(XPoint *) malloc((unsigned) 2*(Ncells_n_1+10)*sizeof(int));
   XSetForeground(display,context,color[196].pixel);
   for(j=0;j<Ncells_n_1;j++){
      (FIBR+j)->x = XtoGaphics((CELL_n_1+j)->x);
      (FIBR+j)->y = DMtoGaphics((CELL_n_1+j)->delta_m);
      XFillArc(display,window_U,context,(FIBR+j)->x-2,(FIBR+j)->y-2,4,4,0,64*360);
   }
   free(FIBR);
   FIBR=(XPoint *) malloc((unsigned) 2*(Ncells_n+10)*sizeof(int));
   for(j=0;j<Ncells_n;j++){
      (FIBR+j)->x = XtoGaphics((CELL_n+j)->x);
      (FIBR+j)->y = DMtoGaphics((CELL_n+j)->delta_m);
      XSetForeground(display,context,
         color[(int)(195*(1.-(((CELL_n+j)->rmax)-ymin)/(ymax-ymin)))].pixel);
      XFillArc(display,window_U,context,(FIBR+j)->x-2,(FIBR+j)->y-2,4,4,0,64*360);
   }
   free(FIBR);

   XSetForeground(display,context,1);
   XDrawRectangle(display,window_U,context,GDx-1,GDy-1,GLx+2,GLy+2);
   for(i=0;i<=10;i++){
      sprintf(Label,"%.2f",(0.1*i*Xmax));
      XDrawString(display,window_U,context,
         XtoGaphics(0.1*i*Xmax)-10,GDy+GLy+16,Label,strlen(Label));
   }
   XSetForeground(display,context,1);
   sprintf(Label,"Chemotactic signaling at cells, delta m");
   XDrawString(display,window_U,context,GDx,GDy-8,Label,strlen(Label));
   XSetForeground(display,context,color[180].pixel);
   sprintf(Label,"PDGF concentration, u");
   XDrawString(display,window_U,context,GDx+GLx-130,GDy-8,Label,strlen(Label));

// *** FIBROBLAST DENSITY, v ********************************************************************************
   XSetForeground(display,context,color[196].pixel);
   XFillRectangle(display,window_V,context,GDx,GDy,GLx+1,GLy+1);

   XSetForeground(display,context,1);
   for(j=0;j<N_V;j++)
       (DENS+j)->y = VtoGaphics(V_n[j]);
   XDrawLines(display,window_V,context,DENS,j,0);
   XDrawRectangle(display,window_V,context,GDx-1,GDy-1,GLx+2,GLy+2);
   sprintf(Label,"Fibroblast density, v");
   XDrawString(display,window_V,context,GDx,GDy-8,Label,strlen(Label));

   XSetForeground(display,context,color[180].pixel);
   for(j=0;j<N_V;j++)
       (DENS+j)->y = RtoGaphics(Average_r_max(j));
   XDrawLines(display,window_V,context,DENS,j,0);
   sprintf(Label,"r max, averaged over cells");
   XDrawString(display,window_V,context,GDx+GLx-130,GDy-8,Label,strlen(Label));

// *** CELL POSITIONS **************************************************************************************
   FIBR=(XPoint *) malloc((unsigned) 2*(Ncells_n+10)*sizeof(int));
   XSetForeground(display,context,color[196].pixel);
   XFillRectangle(display,window_C,context,GDx-3,GDy-3,GLx+6,GLy+6);
   XSetForeground(display,context,1);
   XDrawRectangle(display,window_C,context,GDx-2,GDy-2,GLx+4,GLy+4);
   for(i=0;i<Ncells_n;i++){
      (FIBR+i)->x = XtoGaphics((CELL_n+i)->x);
      (FIBR+i)->y = CtoGaphics((CELL_n+i)->y);
      XSetForeground(display,context,
         color[(int)(195.*(1.-(((CELL_n+i)->rmax)-ymin)/(ymax-ymin)))].pixel);
      XFillArc(display,window_C,context,(FIBR+i)->x-2,(FIBR+i)->y-2,4,4,0,64*360);
      XDrawArc(display,window_C,context,(FIBR+i)->x-2,(FIBR+i)->y-2,4,4,0,64*360);
   }
   free(FIBR);
}

void ColorScale(){
int i;
Arg args[5];
   for(i=0;i<196;i++){
      color[i].pixel=i;
      color[i].blue=32000+32000.0*tanh((i-120)*0.01);
      color[i].green=32000+32000.0*cos(i*4.0*3.1416/250.0);
      color[i].red=32000+32000.0*cos(i*2.0*3.1416/250.0);
      color[i].flags=DoRed|DoGreen|DoBlue;
      XAllocColor(display,DefaultColormap(display,0),color+i);
   }
   color[196].pixel=196;
   color[196].blue=64000;
   color[196].green=64000;
   color[196].red=64000;
   color[196].flags=DoRed|DoGreen|DoBlue;
   XAllocColor(display,DefaultColormap(display,0),color+196);
}

void Print_Time(){
char number[50];
Arg args[5];

   sprintf(number,"time [h]: %.1f",TIME);
   XtSetArg(args[0],XtNlabel,number);
   XtSetValues(iter_box,args,1);
}

void Print_Data(Widget w, XtPointer call_data, XtPointer client_data){
char filename[50];
int i;

   FILE *file_U;
   sprintf(filename,"%s/run_%05d_U",coeff,i_run);
   file_U=fopen(filename,"w");
   fprintf(file_U,"\n %d",N_U);
   for(i=0;i<N_U;i++)
      fprintf(file_U,"\n %f   %f",X_u[i],U_n[i]);
   fclose(file_U);

   FILE *file_V;
   sprintf(filename,"%s/run_%05d_V",coeff,i_run);
   file_V=fopen(filename,"w");
   fprintf(file_V,"\n %d",N_V);
   for(i=0;i<N_V;i++)
      fprintf(file_V,"\n %f   %f",X_v[i],V_n[i]);
   fclose(file_V);

   FILE *file_av;
   sprintf(filename,"%s/run_%05d_rav",coeff,i_run);
   file_av=fopen(filename,"w");
   fprintf(file_av,"\n %d",N_V);
   for(i=0;i<N_V;i++)
      fprintf(file_av,"\n %f  %f  %f  %f  %f  %f  %f",X_v[i],Average_r_max(i),Average_r_i(i),
	Average_e_max(i),Average_e_i(i),Average_p_max(i),Average_Delta_m(i));
   fclose(file_av);

   FILE *file_rmax;
   sprintf(filename,"%s/run_%05d_rmax",coeff,i_run);
   file_rmax=fopen(filename,"w");
   fprintf(file_rmax,"\n %d",Ncells_n);
   for(i=0;i<=Ncells_n;i++)
      fprintf(file_rmax,"\n %f   %f",(CELL_n+i)->x,(CELL_n+i)->rmax);
   fclose(file_rmax);

   PenDepth+=VfrontPosition(0.5)-x_dermis;
   CellsInvaded+=(float)(CellsInTheClot())/Ncells_0;
}

void Average_and_print_depth(){
char filename[50];
FILE *file_Vfront;

   sprintf(filename,"%s/%d.out",directory,(int)(TIME));
   file_Vfront=fopen(filename,"w");
   fprintf(file_Vfront,"\n %f   %f",PenDepth/n_runs,CellsInvaded/n_runs);
   fclose(file_Vfront);
}

void CheckMaxLevels(){
double max,ymax,y;
int i;

   max=U_n[0];
   for(i=1;i<N_U;i++){
      if(max<U_n[i])
         max=U_n[i];
   }
   if(max!=0.){
      y=pow(10.,(int)(0.43429448190325176*log(max)));
      ymax=y*((int)(max/y)+1);
      if(ymax!=Umax){
         Umax=ymax;
         DrawFrame(Rmin,Rmax);
         PlotData(Rmin,Rmax);
      }
   }

   max=(CELL_n+0)->delta_m;
   for(i=1;i<=Ncells_n;i++){
      if(max<(CELL_n+i)->delta_m)
         max=(CELL_n+i)->delta_m;
   }
   if(max!=0.){
      y=pow(10.,(int)(0.43429448190325176*log(max)));
      ymax=y*((int)(max/y)+1);
      if(ymax!=DMmax){
         DMmax=ymax;
         DrawFrame(Rmin,Rmax);
         PlotData(Rmin,Rmax);
      }
   }

   max=V_n[0];
   for(i=1;i<N_V;i++){
      if(max<V_n[i])
         max=V_n[i];
   }
   if(max!=0.){
      y=pow(10.,(int)(0.43429448190325176*log(max)));
      ymax=y*((int)(max/y)+1);
      if(ymax!=Vmax){
         Vmax=ymax;
         DrawFrame(Rmin,Rmax);
         PlotData(Rmin,Rmax);
      }
   }
}

void PrintMovieData(){
FILE *movie;
char filename[50];
int i;

   sprintf(filename,"Movie/%05.1f",TIME);
   movie=fopen(filename,"w");
   for(i=0;i<Ncells_n;i++)
      fprintf(movie,"%f   %f   %f\n",(CELL_n+i)->x,(CELL_n+i)->y,(CELL_n+i)->rmax);
   fclose(movie);
}