#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <X11/StringDefs.h>
#include <X11/Xaw/Cardinals.h>
#include <X11/Xaw/Command.h>
#include <X11/Xaw/Box.h>
#include <X11/Xaw/List.h>
#include <assert.h>

#define pi				3.14159265358979
#define sqrtPI          sqrt(pi)
#define DIMENSION		600
#define MIN(x,y)		((x<y) ? x:y)
#define nu				0.01

Widget toplevel,time_box,size_box,alpha_box;
GC cont;
Window win;
Display *disp;
XColor color[10];
static XtAppContext app_con;

#define GRID_POINTS    20000	// Required for evaluation of Error function
typedef struct {double XMIN, XMAX, DX; double Table[GRID_POINTS][2];} LookUpTable;
LookUpTable ERFCX, *ptr_erfcx=&ERFCX, ERF, *ptr_erf=&ERF;
double EVALUATE(double,LookUpTable *POINT,int);
void GET_DATA(LookUpTable *erf, LookUpTable *erfcx);

double GAUSS(),RAN(),Rotation_angle(double,double),PDF_React(double,double,double),
		PDF_Reflect(double,double),P_first(double,double);
void Initialize_receptors(),PeriodicBC(),Read_parameters(char []),Initialize_Ras(),
		Initialize_receptors(),Update_receptors(double),Count_states(),
		Gauss_propagation(double,double,double,double,double,double *,double *,double *,double *),
		PeriodicBC(int),Circle_Point_Picking(double,double *,double *),Update_average_values(),
		Ras_distance_to_closest_receptor(int),Ras_reaction_diffusion();

void Quit(Widget,XtPointer,XtPointer),Start(Widget,XtPointer,XtPointer),Draw(),ColorScale(),
	Advance_Step(Widget,XtPointer,XtPointer),Switch_Graphics(Widget,XtPointer,XtPointer);

struct particles {
double x;
double y;
short state;	// active/inactive state
double dmin;	// distance to the closest receptor
int n;			// index of the closest receptor
};
struct receptors {
double x;
double y;
short state;	// enzyme-on/off state
};
struct particles *Ras, *B_Ras;
struct receptors *Rec;

int N_Rec,N_Rec_enz,N_Rec_off,N_Ras;
double Lxy,eta_RE,S,SIGMA,Diff,K_RE,TIME,k_r_on,k_r_off,t_last,t_out,N_act,dS,TSTEP,k_in,
	kappa,Da,delta_t,alpha,N_in_sum,N_RE_sum,n_steps;
int on_off=0,show=1,enzyme_on_off;
FILE *OUTPUT;
char NAMEIN[100],NAMEOUT[100];

// *****************************************************************************
// ******* INITIALIZATION ***** (Function main() is moved to the end) **********
// *****************************************************************************
void Initialize_receptors(){
int i,n,overlaps;
double x,y,dx,dy;

   Lxy=S*sqrt(N_Rec/eta_RE);
   N_Rec_enz=0;
   Rec=(struct receptors *) malloc( 2*N_Rec*(2*sizeof(double)+sizeof(int)) );

   for(i=0;i<N_Rec;i++){
      if(enzyme_on_off==1)
         Rec[i].state = 0;
	  else
         Rec[i].state = 1;
      N_Rec_enz+=Rec[i].state;
   }
   N_Rec_off=N_Rec-N_Rec_enz;

   if(N_Rec>1){
      n=0;
      do{
         x=Lxy*(0.05+0.9*RAN());
         y=Lxy*(0.05+0.9*RAN());
	     overlaps=0;
	     for(i=n-1;i>=0;i--){
		    dx=Rec[i].x-x;				// evaluate a vector of distances
		    dy=Rec[i].y-y;
	        if( sqrt(dx*dx+dy*dy)<2*S )	// avoid overlaps with other particles
		       overlaps=1;
	     }
         if(overlaps==0){
            Rec[n].x=x;
            Rec[n].y=y;
		    n++;
		 }
      }while(n<N_Rec);
   }else{
      Rec[0].x=0.5*Lxy;
      Rec[0].y=0.5*Lxy;
   }
}
void Initialize_Ras(){
int i,j,m;
double dmin,dist;

   Ras=(struct particles *) malloc(2*N_Ras*(3*sizeof(double)+sizeof(int)+sizeof(short)));
   
   for(i=0;i<N_Ras;i++){
      Ras[i].state=0;
      do{
	     Ras[i].x=Lxy*RAN();
	     Ras[i].y=Lxy*RAN();

	     dmin=1.e10;
         for(j=0;j<N_Rec;j++){
            dist=sqrt((Ras[i].x-Rec[j].x)*(Ras[i].x-Rec[j].x)+
			          (Ras[i].y-Rec[j].y)*(Ras[i].y-Rec[j].y))-S;
	        if(dmin>dist){
	           dmin=dist;
	           m=j;
	        }
	     }
      }while(dmin < 0.);
      Ras[i].dmin=dmin;
      Ras[i].n=m;
   }
}
void Read_parameters(char inputfile[20]){
FILE *inputpar;
char sym[20];

   inputpar=fopen(inputfile,"r");
   fscanf(inputpar,"%s %lf",	&sym,&eta_RE);
   fscanf(inputpar,"%s %d",		&sym,&N_Rec);
   fscanf(inputpar,"%s %lf",	&sym,&S);
   fscanf(inputpar,"%s %d",		&sym,&N_Ras);
   fscanf(inputpar,"%s %d",		&sym,&enzyme_on_off);
   fscanf(inputpar,"%s %lf",	&sym,&k_r_on);
   fscanf(inputpar,"%s %lf",	&sym,&k_r_off);
   fscanf(inputpar,"%s %lf",	&sym,&Da);
   fscanf(inputpar,"%s %lf",	&sym,&Diff);
   // Diff=2*10^5 nm^2/s --> 2*10^-9 cm^2/s
   // Diff=2*10^4 nm^2/s --> 2*10^-10 cm^2/s
   // Diff=2*10^3 nm^2/s --> 2*10^-11 cm^2/s
   fscanf(inputpar,"%s %lf",	&sym,&kappa);
   fscanf(inputpar,"%s %lf",	&sym,&TSTEP);
   fscanf(inputpar,"%s %lf",	&sym,&dS);
   fscanf(inputpar,"%s %lf",	&sym,&t_out);

   SIGMA   = sqrt(4.*Diff*TSTEP);
   K_RE=kappa*Diff/(2.*pi*S);
   k_in=Da*Diff/(S*S);
   printf("\n k_in=%e,  K_RE=%e\n",k_in,K_RE);
   fclose(inputpar);
}

// *****************************************************************************
// ******* BASIC FUNCTIONS OF THE ALGORITHM ************************************
// *****************************************************************************
Boolean Main_Step(XtPointer client_data){
int i;
char str[100];
Arg args[5];
double N_RE_av,N_in_av;

   switch(on_off){
      case 1:
      // *** Update receptors ***
      Ras_reaction_diffusion();
      Update_average_values();
      TIME+=delta_t;

      // *** Output ***
      if(TIME-t_last>t_out){
         if(show==1)
	        Draw();
	     i=0;
         XtSetArg(args[i],XtNwidth,200); i++;
         sprintf(str,"Time: %.2e sec",TIME);
         XtSetArg(args[i],XtNlabel,str); i++;
         XtSetValues(time_box,args,i);

         N_in_av=N_in_sum/n_steps;
         N_RE_av=N_RE_sum/n_steps;
	     if(N_in_av==0. || N_RE_av==0.)
	        alpha=0.;
         else
	        alpha=N_act*Lxy*Lxy/(N_in_av*N_RE_av*TIME*Diff);
	     i=0;
         XtSetArg(args[i],XtNwidth,150); i++;
         sprintf(str,"alpha=%.2f",alpha);
         XtSetArg(args[i],XtNlabel,str); i++;
         XtSetValues(alpha_box,args,i);
	  
	     Count_states();
	     fprintf(OUTPUT,"%f  %d  %d  %f  %f  %f\n",
				TIME,N_Rec_enz,N_Rec_off,N_RE_av,N_in_av,alpha);
	     t_last=TIME;
      }
      return(False);
      break;
      case 0:
      return(True);
      break;
   }
}

void Update_average_values(){
int i,n;

   n=0;
   for(i=0;i<N_Rec;i++)
      n+=Rec[i].state;	// Count enzyme-bound receptors
   N_RE_sum+=1.*n;
   n=0;
   for(i=0;i<N_Ras;i++)
      n+=Ras[i].state;	// Count active Ras
   N_in_sum+=1.*(N_Ras-n);

   n_steps+=1.;
}
void Count_states(){
int i;

   N_Rec_enz=0;
   for(i=0;i<N_Rec;i++){
      if(Rec[i].state==1)
         N_Rec_enz++;
   }
   N_Rec_off=N_Rec-N_Rec_enz;
}

void Ras_reaction_diffusion(){
double p_x,p_y,dr,dx,dy,z,z0,P,dt[5];
int i,n,n_react;

   // *** Find the minimum Ras-receptor distance ***
   dr=1.e10;
   n=0;
   for(i=0;i<N_Ras;i++){
      PeriodicBC(i);
	  Ras_distance_to_closest_receptor(i);
	  if(dr>Ras[i].dmin)
	     dr=Ras[i].dmin;
	  n+=Ras[i].state;
   }
   // *** Find minimum diffusion timestep
   if(dr<dS)
      dt[0]=TSTEP;
   else
      dt[0]=dr*dr/(4.*Diff);

   // IMPORTANT: the following reactions happen independently within dlta_t!
   n_react=1;
   // *** Timestep for sampling of inactivation reaction ***
   if(n>0 && k_in>0.){
      dt[n_react]=nu/k_in;
      n_react++;
   }
   if(enzyme_on_off==1){
      // *** Timestep for sampling of enzyme binding ***
      n=0;
      for(i=0;i<N_Rec;i++)
         n+=Rec[i].state;
      if((N_Rec-n)>0 && k_r_on>0.){
         dt[n_react]=nu/k_r_on;
         n_react++;
      }
      // *** Timestep for sampling of enzyme dissociation ***
      if(n>0 && k_r_off>0.){
         dt[n_react]=nu/k_r_off;
         n_react++;
      }
   }
   // *** Find min. time step for the whole reactive system
   delta_t=1.0e10;
   for(i=0;i<n_react;i++){
      if(delta_t>dt[i])
         delta_t=dt[i];
   }

   if(enzyme_on_off==1)
      Update_receptors(delta_t);

   for(i=0;i<N_Ras;i++){
      if(Ras[i].dmin<dS){
         n=Ras[i].n;
         p_x=Ras[i].x;
		 p_y=Ras[i].y;
		 Gauss_propagation(Rec[n].x,Rec[n].y,p_x,p_y,delta_t,&Ras[i].x,&Ras[i].y,&z0,&z);
		 if( Ras[i].state==0 ){
            P=1.-PDF_React(z,z0,K_RE)/PDF_Reflect(z,z0);
			if( Rec[n].state==1 && P>RAN() ){
			   Ras[i].state=1;
			   N_act+=1.;
			}
		 }else{
		    if( P_first(k_in,delta_t)>RAN() )
			   Ras[i].state=0;
		 }
      }else{
         if( Ras[i].state==1 && P_first(k_in,delta_t)>RAN() )
            Ras[i].state=0;
         dr=sqrt(delta_t*4.*Diff);
         // *** Update coordinates ***
		 Circle_Point_Picking(dr,&dx,&dy);
		 Ras[i].x += dx;
		 Ras[i].y += dy;
      }
   }
}
void Update_receptors(double dt){
int i;

   for(i=0;i<N_Rec;i++){
      if(Rec[i].state==0){
         if( P_first(k_r_on,dt)>RAN() )
		    Rec[i].state=1;
	  }else{
         if( P_first(k_r_off,dt)>RAN() )
		    Rec[i].state=0;
      }
   }
}
void Ras_distance_to_closest_receptor(int j){
int i,n;
double d,dmin,dx,dy;

   dmin=1.e10;
   for(i=0;i<N_Rec;i++){
      d=sqrt((Ras[j].x-Rec[i].x)*(Ras[j].x-Rec[i].x)+(Ras[j].y-Rec[i].y)*(Ras[j].y-Rec[i].y))-S;
      if(dmin>d){
         dmin=d;
         n=i;
      }
   }
   if(dmin<0.){
	  dmin=RAN()*dS;
	  Circle_Point_Picking(S+dmin,&dx,&dy);
	  Ras[j].x=Rec[n].x+dx;
	  Ras[j].y=Rec[n].y+dy;
   }
   Ras[j].dmin=dmin;
   Ras[j].n=n;
}

void Circle_Point_Picking(double r, double *xr, double *yr){
double x1,x2,x12,x22,sum;

   do{
      x1=(1.-2.*RAN());   x12=x1*x1;
      x2=(1.-2.*RAN());   x22=x2*x2;
      sum=x12+x22;
   }while(sum>=1.);
   *xr=r*(x12-x22)/sum;
   *yr=r*2.*x1*x2/sum;
}
void Gauss_propagation(double Xr, double Yr, double Xi, double Yi, double DT,
		double *Xnew, double *Ynew, double *DR_curr, double *DR_new){
double ksi,eta,rho,alpha,beta,Ra;

	Ra=sqrt((Xi-Xr)*(Xi-Xr)+(Yi-Yr)*(Yi-Yr));
	*DR_curr = Ra - S;	
	ksi=fabs( *DR_curr+GAUSS()*sqrt(2.*DT*Diff) );
	eta=GAUSS()*sqrt(2.*DT*Diff);
	rho=sqrt(ksi*ksi+eta*eta);
	alpha=Rotation_angle(ksi,eta);
	beta=Rotation_angle((Xi-Xr),(Yi-Yr));
	*Xnew = Xr + S*cos(beta) + rho*cos(alpha+beta);
	*Ynew = Yr + S*sin(beta) + rho*sin(alpha+beta);
	*DR_new=sqrt((*Xnew-Xr)*(*Xnew-Xr)+(*Ynew-Yr)*(*Ynew-Yr)) - S;
}
double Rotation_angle(double DX, double DY){
double angle;
   if(fabs(DX)<1.e-10){
      if(DY>0.)
         angle=pi/2.;
      else
         angle=pi+pi/2.;
   }
   else{
      if(DX>0. && DY>0.)
         angle=atan(DY/DX);
      if(DX<0. && DY>0.)
         angle=pi-atan(DY/(-DX));
      if(DX<0. && DY<0.)
         angle=pi+atan(DY/DX);
      if(DX>0. && DY<0.)
         angle=2.*pi-atan(-DY/DX);
   }
   return angle;
}
void PeriodicBC(int J){
double X,Y;

   X=Ras[J].x;
   Y=Ras[J].y;
   if(X>0. && X<Lxy){
   }else{
      if( X>Lxy )
	     Ras[J].x = (double) fmod(X,Lxy);
	  else
		 Ras[J].x = Lxy + (double) fmod(X,Lxy);
   }
   if(Y>0. && Y<Lxy){
   }else{ 
      if( Y>Lxy )
	     Ras[J].y = (double) fmod(Y,Lxy);
	  else
		 Ras[J].y = Lxy + (double) fmod(Y,Lxy);
   }
}

double P_first(double k, double dt){
   return 1.-exp(-k*dt);
}
double PDF_React(double x, double x0, double k){
double z,z0,f1,f2,mu_plus;
   z  = x/SIGMA;
   z0 = x0/SIGMA;
   mu_plus = 2.*k*TSTEP/sqrt(4.*Diff*TSTEP);
   f1 = ( exp(-(z-z0)*(z-z0)) + exp(-(z+z0)*(z+z0)) )/(SIGMA*sqrtPI); 
   f2 = (k/Diff)*exp(-(z+z0)*(z+z0))*EVALUATE((z+z0+mu_plus),ptr_erfcx,0);
   return f1 - f2;
}
double PDF_Reflect(double x, double x0){
double z,z0;
   z  = x/SIGMA;
   z0 = x0/SIGMA;
   return ( exp(-(z-z0)*(z-z0)) + exp(-(z+z0)*(z+z0)) )/(SIGMA*sqrtPI); 
}

double RAN(){
   return drand48();
}
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;
   }
}

void GET_DATA(LookUpTable *erf, LookUpTable *erfcx){
int I, index;
FILE *ifp; 
double a, b, c;    

   assert( (ifp=fopen("ERFCX.dat","r")) != NULL);
   fscanf(ifp,"%lf%lf", &(erfcx -> XMIN), &(erfcx -> XMAX));
   fscanf(ifp,"%lf%lf", &(erfcx -> DX),   &(erfcx -> DX));
   a = (erfcx -> XMAX); 
   b = (erfcx -> XMIN);
   c = (erfcx -> DX);
   I = (int) ((a-b)/c);
   for(index=0; index<=I; index++)
	  fscanf(ifp, "%lf%lf", &((*erfcx).Table)[index][0], &((*erfcx).Table)[index][1]);
   fclose(ifp);

   assert( (ifp=fopen("ERF.dat","r")) != NULL);
   fscanf(ifp,"%lf%lf", &(erf -> XMIN), &(erf -> XMAX));
   fscanf(ifp,"%lf%lf", &(erf -> DX),   &(erf -> DX));
   a = (erf -> XMAX); 
   b = (erf -> XMIN);
   c = (erf -> DX);
   I = (int) ((a-b)/c);
   for(index=0; index<=I; index++)
	  fscanf(ifp, "%lf%lf", &((*erf).Table)[index][0], &((*erf).Table)[index][1]);
   fclose(ifp);
   return;
}
double EVALUATE(double x, LookUpTable *POINT, int MODE){
double k, b;
int IND;

   assert( x>=(*POINT).XMIN );
   if( x<=(*POINT).XMAX ){ 
	  IND = (int) (floor(x/(*POINT).DX));
      k = (POINT -> Table)[IND][0];
      b = (POINT -> Table)[IND][1];
      return (k*x + b); 
   }else{
      if(MODE)			/* asymptotic expressions*/
         return 1.0;	/* error function*/
	  else 
         return 1/(sqrtPI*x);  /* scaled complementary error function*/
   }
}


// *****************************************************************************
// ******* MAIN() **************************************************************
// *****************************************************************************
main(int argc,char *argv[]){
Widget image_box,control_box,quit,start,main_box,step,draw;
Arg args[10];
int i;
static String calc_on_off[]={"STOP","GO",NULL};
static String display_on_off[]={"Graphics Off","Graphics On",NULL};

   sprintf(NAMEIN,"%s",argv[1]);
   sprintf(NAMEOUT,"%s",argv[2]);
   OUTPUT=fopen(NAMEOUT,"w");
   GET_DATA(ptr_erf, ptr_erfcx);

   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++;
   XtSetArg(args[i],XtNwidth,50+DIMENSION); i++;
   XtSetArg(args[i],XtNheight,50+DIMENSION); i++;
   main_box=XtCreateManagedWidget("main_box", boxWidgetClass,toplevel,args,i);
   i=0;
   XtSetArg(args[i],XtNwidth,DIMENSION); i++;
   control_box=XtCreateManagedWidget("control_box", boxWidgetClass,main_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlabel,"Initialize data"); i++;
   start=XtCreateManagedWidget("start",commandWidgetClass,control_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlist,calc_on_off); i++;
   XtSetArg(args[i],XtNwidth,100); i++;
   step=XtCreateManagedWidget("step",listWidgetClass,control_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlist,display_on_off); i++;
   XtSetArg(args[i],XtNwidth,200); i++;
   draw=XtCreateManagedWidget("draw",listWidgetClass,control_box,args,i);
   i=0;
   XtSetArg(args[i],XtNlabel,"Quit"); i++;
   quit=XtCreateManagedWidget("quit",commandWidgetClass,control_box,args,i);
   i=0;
   XtSetArg(args[i],XtNwidth,DIMENSION); i++;
   XtSetArg(args[i],XtNheight,DIMENSION); i++;
   image_box=XtCreateManagedWidget("image_box",simpleWidgetClass,main_box,args,i);
   i=0;
   XtSetArg(args[i],XtNwidth,200); i++;
   XtSetArg(args[i],XtNlabel,"Box size:"); i++;
   size_box=XtCreateManagedWidget("Box size",labelWidgetClass,main_box,args,i);
   i=0;
   XtSetArg(args[i],XtNwidth,200); i++;
   XtSetArg(args[i],XtNlabel,"Time:"); i++;
   time_box=XtCreateManagedWidget("time",labelWidgetClass,main_box,args,i);
   i=0;
   XtSetArg(args[i],XtNwidth,150); i++;
   XtSetArg(args[i],XtNlabel,"alpha=0"); i++;
   alpha_box=XtCreateManagedWidget("alpha",labelWidgetClass,main_box,args,i);

   XtAddCallback(start,XtNcallback,Start,NULL);
   XtAddCallback(step,XtNcallback,Advance_Step,NULL);
   XtAddCallback(draw,XtNcallback,Switch_Graphics,NULL);
   XtAddCallback(quit,XtNcallback,Quit,NULL);

   XtRealizeWidget(toplevel);
   disp=XtDisplay(image_box);
   win=XtWindow(image_box);
   cont=XCreateGC(disp,win,0,NULL);
   ColorScale();
   XtAppMainLoop(app_con);
}

// *****************************************************************************
// ******* G U I ***************************************************************
// *****************************************************************************
void Start(Widget w, XtPointer client_data, XtPointer call_data){
char str[50];
Arg args[10];
int i;

   TIME=0.;
   t_last=0.;
   N_act=0.;
   N_in_sum=0.;
   N_RE_sum=0.;
   n_steps=0.;
   Read_parameters(NAMEIN);
   Initialize_receptors();
   Initialize_Ras();
   Draw();
   i=0;
   XtSetArg(args[i],XtNwidth,200); i++;
   sprintf(str,"Box size: %.2f nm",Lxy);
   XtSetArg(args[i],XtNlabel,str); i++;
   XtSetValues(size_box,args,i);
   i=0;
   XtSetArg(args[i],XtNwidth,200); i++;
   sprintf(str,"Time: %.2e sec",TIME);
   XtSetArg(args[i],XtNlabel,str); i++;
   XtSetValues(time_box,args,i);

   fclose(OUTPUT);
   OUTPUT=fopen(NAMEOUT,"w");

   XtAppAddWorkProc(app_con,Main_Step,NULL);
}
void Quit(Widget w, XtPointer call_data, XtPointer client_data){
   fclose(OUTPUT);
   exit(0);
}
void Advance_Step(Widget w, XtPointer client_data, XtPointer call_data){
   XawListReturnStruct *item = (XawListReturnStruct*)call_data;
   on_off=item->list_index;
   XtAppAddWorkProc(app_con,Main_Step,NULL);
}
void Switch_Graphics(Widget w, XtPointer client_data, XtPointer call_data){
   Draw();
   XawListReturnStruct *item = (XawListReturnStruct*)call_data;
   show=item->list_index;
}

// *****************************************************************************
// ******* GRAPHICS ************************************************************
// *****************************************************************************
void Draw(){
int i,I,J,D_Rec,D_Ras=4;
char Label[10];

   XClearWindow(disp,win);
   XSetForeground(disp,cont,color[0].pixel);
   XFillRectangle(disp,win,cont,0,0,DIMENSION,DIMENSION);

   for(i=0;i<N_Rec;i++){
      D_Rec=(int)(DIMENSION*2*S/Lxy+0.5);
	  if(D_Rec<10)
	     D_Rec=10;

      I=(int)(DIMENSION*Rec[i].x/Lxy-0.5*D_Rec+0.5);
      J=(int)(DIMENSION*Rec[i].y/Lxy-0.5*D_Rec+0.5);
	  if(Rec[i].state==0)
         XSetForeground(disp,cont,color[1].pixel);
	  else
         XSetForeground(disp,cont,color[4].pixel);
      XFillArc(disp,win,cont,I,J,D_Rec,D_Rec,0,64*360);
      XSetForeground(disp,cont,1);
      XDrawArc(disp,win,cont,I,J,D_Rec,D_Rec,0,64*360);
   }
   for(i=0;i<N_Ras;i++){
      I=(int)(DIMENSION*Ras[i].x/Lxy-0.5*D_Ras+0.5);
      J=(int)(DIMENSION*Ras[i].y/Lxy-0.5*D_Ras+0.5);
	  if(Ras[i].state==0){
         XSetForeground(disp,cont,color[3].pixel);
         XFillArc(disp,win,cont,I,J,D_Ras,D_Ras,0,64*360);
         XDrawArc(disp,win,cont,I,J,D_Ras,D_Ras,0,64*360);
	  }else{
         XSetForeground(disp,cont,color[2].pixel);
         XFillArc(disp,win,cont,I,J,D_Ras,D_Ras,0,64*360);
         XDrawArc(disp,win,cont,I,J,D_Ras,D_Ras,0,64*360);
	  }
   }
}
void ColorScale(){
int i;
// *** WHITE ***
   i=0;
   color[i].pixel=i;	color[i].blue=64000;	color[i].green=64000;	color[i].red=64000;
   color[i].flags=DoRed|DoGreen|DoBlue;
   XAllocColor(disp,DefaultColormap(disp,0),color+i);
// *** YELLOW ***
   i=1;
   color[i].pixel=i;	color[i].blue=0;	color[i].green=64000;	color[i].red=64000;
   color[i].flags=DoRed|DoGreen|DoBlue;
   XAllocColor(disp,DefaultColormap(disp,0),color+i);
// *** RED ***
   i=2;
   color[i].pixel=i;	color[i].blue=0;	color[i].green=0;	color[i].red=64000;
   color[i].flags=DoRed|DoGreen|DoBlue;
   XAllocColor(disp,DefaultColormap(disp,0),color+i);
// *** GREEN ***
   i=3;
   color[i].pixel=i;	color[i].blue=0;	color[i].green=64000;	color[i].red=0;
   color[i].flags=DoRed|DoGreen|DoBlue;
   XAllocColor(disp,DefaultColormap(disp,0),color+i);
// *** BLUE ***
   i=4;
   color[i].pixel=i;	color[i].blue=64000;	color[i].green=0;	color[i].red=0;
   color[i].flags=DoRed|DoGreen|DoBlue;
   XAllocColor(disp,DefaultColormap(disp,0),color+i);
}

// *****************************************************************************
// ******* THE END *************************************************************
// *****************************************************************************