}
     for(k=(length+1); k<=MAXL_N; k++)
     {
      temp_seq[k]=0;
     }
     length=sq_length[xd5][yd5];
     sq_length[xd4][yd4]=sq_length[xd5][yd5];
     t_n[xd4][yd4]=t_n[xd5][yd5]; r_n[xd4][yd4]=r_n[xd5][yd5]; q_n[xd4][yd4]=q_n[xd5][yd5];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd4][yd4][k]=sq_seq[xd5][yd5][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd4][yd4][k]=0;
     }
     length=temp_l;
     sq_length[xd5][yd5]=temp_l;
     t_n[xd5][yd5]=temp_t; r_n[xd5][yd5]=temp_r; q_n[xd5][yd5]=temp_q;
     for(k=1; k<=length; k++)
     {
      sq_seq[xd5][yd5][k]=temp_seq[k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd5][yd5][k]=0;
     }
     length=sq_length[xd6][yd6];
     temp_l=length;
     temp_t=t_n[xd6][yd6]; temp_r=r_n[xd6][yd6]; temp_q=q_n[xd6][yd6];
     for(k=1; k<=length; k++)
     {
      temp_seq[k]=sq_seq[xd6][yd6][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      temp_seq[k]=0;
     }
     length=sq_length[xd1][yd1];
     sq_length[xd6][yd6]=sq_length[xd1][yd1];
     t_n[xd6][yd6]=t_n[xd1][yd1]; r_n[xd6][yd6]=r_n[xd1][yd1]; q_n[xd6][yd6]=q_n[xd1][yd1];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd6][yd6][k]=sq_seq[xd1][yd1][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd6][yd6][k]=0;
     }
     length=sq_length[xd3][yd3];
     sq_length[xd1][yd1]=sq_length[xd3][yd3];
     t_n[xd1][yd1]=t_n[xd3][yd3]; r_n[xd1][yd1]=r_n[xd3][yd3]; q_n[xd1][yd1]=q_n[xd3][yd3];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd1][yd1][k]=sq_seq[xd3][yd3][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd1][yd1][k]=0;
     }
     length=temp_l;
     sq_length[xd3][yd3]=temp_l;
     t_n[xd3][yd3]=temp_t; r_n[xd3][yd3]=temp_r; q_n[xd3][yd3]=temp_q;
     for(k=1; k<=length; k++)
     {
      sq_seq[xd3][yd3][k]=temp_seq[k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd3][yd3][k]=0;
     }
  }
 
 void mc_step(void)
  {
   mc_stepi();

   for (j=1; j<=SIZE; j++)
   {
    xy_choose();
    if(randd()<DIFFRATE)
    {
     toff_marg();
    }

    length=sq_length[x][y];

    if(!length) {continue;} /* if empty site -> continue */

    if(randd()<DEATH)        /* if death occurs -> overwrite, continue */
    {
     for(k=1; k<=length; k++)     {      sq_seq[x][y][k]=0;     }     t_n[x][y]=0;     r_n[x][y]=0;     q_n[x][y]=0;     sq_length[x][y]=0;
     n_denz=(n_denz*SIZE-length)/SIZE;     continue;
    }

    if (length<CATCH) {continue;}  /* if too short to be a template -> continue */

    neighbourhood();

    if (!(life*empty)) {continue;} /* if there is no at least one      */
				                   /* empty and one filled site        */
				                   /* among the neighbours -> continue */
    rd=0.3333*life;
    if(randd()>rd) {continue;}

    s1=randl(life)+1;
    s1=lifes[s1];
    s2=randl(empty)+1;
    s2=empties[s2];

    if (s1==1) { xs1=left[x];  ys1=y;       }
    if (s1==2) { xs1=right[x]; ys1=y;       }
    if (s1==3) { xs1=x;        ys1=up[y];   } /* coordinates of replicase */
    if (s1==4) { xs1=x;        ys1=down[y]; }

    if (s2==1) { xs2=left[x];  ys2=y;       }
    if (s2==2) { xs2=right[x]; ys2=y;       }
    if (s2==3) { xs2=x;        ys2=up[y];   } /* coordinates of new chain    */
    if (s2==4) { xs2=x;        ys2=down[y]; }


    /* longest continuous sequences of monomers in replicase */

    max_B=r_n[xs1][ys1];
    max_C=q_n[xs1][ys1];

    r=fgv_replicase[max_B];    q=fgv_fidelity[max_C];

    /* longest continuous sequences of monomers in template */

    max_A=t_n[x][y];

    t=fgv_template[max_A];    length;

    rd=randd();

    n_factor=(NUKLPAR-n_denz)/NUKLPAR;    if(rd<n_factor*(t*r/length))                       /* replication */
    {
     sq_length[xs2][ys2]=length;
     for(k=1; k<=length; k++)
     {
      rd=randd();
      if(rd<q)
      {
       sq_seq[xs2][ys2][k]=sq_seq[x][y][k];
      } else
      {
       nukl=sq_seq[x][y][k];       rl=randl(3)+1;       switch (nukl)       {       case 1:        if(rl==1) sq_seq[xs2][ys2][k]=2;        if(rl==2) sq_seq[xs2][ys2][k]=3;        if(rl==3) sq_seq[xs2][ys2][k]=4;        break;       case 2:        if(rl==1) sq_seq[xs2][ys2][k]=1;        if(rl==2) sq_seq[xs2][ys2][k]=3;        if(rl==3) sq_seq[xs2][ys2][k]=4;       break;       case 3:        if(rl==1) sq_seq[xs2][ys2][k]=1;        if(rl==2) sq_seq[xs2][ys2][k]=2;        if(rl==3) sq_seq[xs2][ys2][k]=4;       break;       case 4:        if(rl==1) sq_seq[xs2][ys2][k]=1;        if(rl==2) sq_seq[xs2][ys2][k]=2;        if(rl==3) sq_seq[xs2][ys2][k]=3;       break;       }      }
     }
     rd=randd();     if(rd<b)     {      length_o=length;      length=randl(length-1)+1;
      sq_length[xs2][ys2]=length;
      if(randd()<0.5)      {       for(k=(length+1); k<=length_o; k++)       {        sq_seq[xs2][ys2][k]=0;
       }      } else      {       for(k=1; k<=length; k++)       {        sq_seq[xs2][ys2][k]=sq_seq[xs2][ys2][length_o-length+k];
       }       for(k=(length+1); k<=length_o; k++)       {        sq_seq[xs2][ys2][k]=0;
       }      }     }     if((rd<(a+b)) && (rd>b))
     {
      length++;
      if(randd()<0.5)      {      sq_length[xs2][ys2]++;
      sq_seq[xs2][ys2][length]=randl(4)+1;
      } else      {      for(k=length; k>=2; k--)      sq_seq[xs2][ys2][k]=sq_seq[xs2][ys2][k-1];
      sq_seq[xs2][ys2][1]=randl(4)+1;      sq_length[xs2][ys2]++;
      }     }
     n_denz=(n_denz*SIZE+length)/SIZE;     trq_update();    }
   }
  }

 void save(void)
 {
  tf2=fopen("pos.dat","w+t");  fprintf(tf2,"%d",position);  fflush(tf2);  fclose(tf2);  position++; }

 void result1(void) {  double elokszama=0;  long int lepesszam;  double a_lanchossz=0, S_lanchossz=0, max_lanchossz=0, min_lanchossz=1000;  double A_atlag=0, B_atlag=0, C_atlag=0, D_atlag=0,         As_atlag=0, Bs_atlag=0, Cs_atlag=0, Ds_atlag=0,         t_atlag=0, r_atlag=0, q_atlag=0, t_r_l_atlag;   double A_szoras=0, B_szoras=0, C_szoras=0, D_szoras=0,         As_szoras=0, Bs_szoras=0, Cs_szoras=0, Ds_szoras=0,         t_szoras=0, r_szoras=0, q_szoras=0;  double A_min=1000, B_min=1000, C_min=1000, D_min=1000,         As_min=1000, Bs_min=1000, Cs_min=1000, Ds_min=1000,         t_min=1000, r_min=1000, q_min=1000;    double A_max=0, B_max=0, C_max=0, D_max=0,         As_max=0, Bs_max=0, Cs_max=0, Ds_max=0,         t_max=0, r_max=0, q_max=0, t_r_l_max=0;  double denzitas=0;  long int dat, len, max, maxA, maxB, maxC, maxD, minA, minB, minC, minD;  long int hossz;  double ar;  FILE * tf;
  FILE * tfd;  lepesszam=i;  for(j=1; j<=SIDE; j++)  {   for(k=1; k<=SIDE; k++)   {    hossz=sq_length[j][k];    if(hossz)     {      denzitas+=1; elokszama+=1;      a_lanchossz+=hossz;     maxA=0, maxB=0, maxC=0, maxD=0;      for(l=1; l<=hossz; l++)     {      if(sq_seq[j][k][l]==1) {A_atlag+=1; maxA++;}      if(sq_seq[j][k][l]==2) {B_atlag+=1; maxB++;}      if(sq_seq[j][k][l]==3) {C_atlag+=1; maxC++;}      if(sq_seq[j][k][l]==4) {D_atlag+=1; maxD++;}     }     max=t_n[j][k];     As_atlag+=max;     if(hossz>=CATCH) t_atlag+=fgv_template[max];     max=r_n[j][k];     Bs_atlag+=max;     if(hossz>=CATCH) r_atlag+=fgv_replicase[max];     max=q_n[j][k];     Cs_atlag+=max;     q_atlag+=fgv_fidelity[max];     len=0; max=0;     for(l=1; l<=hossz; l++)     {      dat=sq_seq[j][k][l];      if(dat==4)      {       len++;      }      if(dat!=4)      {       if(len>max) max=len;       len=0;      }     }     Ds_atlag+=max;    }   }  }  a_lanchossz=a_lanchossz/elokszama;  A_atlag=A_atlag/elokszama;  B_atlag=B_atlag/elokszama;  C_atlag=C_atlag/elokszama;  D_atlag=D_atlag/elokszama;  As_atlag=As_atlag/elokszama;  Bs_atlag=Bs_atlag/elokszama;  Cs_atlag=Cs_atlag/elokszama;  Ds_atlag=Ds_atlag/elokszama;  t_atlag=t_atlag/elokszama;  r_atlag=r_atlag/elokszama;  q_atlag=q_atlag/elokszama;  denzitas=denzitas/SIZE;  tf=fopen("result1.dat","a+t");   fprintf(tf,"%10d", lepesszam);   fprintf(tf,"%8.2f",a_lanchossz);   fprintf(tf,"%6.2f",A_atlag);   fprintf(tf,"%6.2f",B_atlag);   fprintf(tf,"%6.2f",C_atlag);   fprintf(tf,"%6.2f",D_atlag);   fprintf(tf,"%6.2f",As_atlag);   fprintf(tf,"%6.2f",Bs_atlag);   fprintf(tf,"%6.2f",Cs_atlag);   fprintf(tf,"%6.2f",Ds_atlag);   fprintf(tf,"%8.4f",t_atlag);   fprintf(tf,"%8.4f",r_atlag);   fprintf(tf,"%8.4f",q_atlag);   fprintf(tf,"%7.3f\n",denzitas);   fflush(tf);  fclose(tf); } void result2(void) {  double elokszama=0;  long int lepesszam;  long int dat, len, max;  restomb1 arany;  restomb2 composition;  long int hossz;  double ar;  FILE * tf;
  FILE * tfd;  lepesszam=i;  for(j=0; j<=MAXL_N; j++) arany[j]=0;  for(j=0; j<=MAXL_N; j++)  {    for(k=0; k<=8; k++) composition[j][k]=0;  }  for(j=1; j<=SIDE; j++)  {   for(k=1; k<=SIDE; k++)   {    hossz=sq_length[j][k];    if(hossz)     {     arany[hossz]+=1;     for(l=1; l<=hossz; l++)     {      if(sq_seq[j][k][l]==1) composition[hossz][1]+=1;      if(sq_seq[j][k][l]==2) composition[hossz][2]+=1;      if(sq_seq[j][k][l]==3) composition[hossz][3]+=1;      if(sq_seq[j][k][l]==4) composition[hossz][4]+=1;     }     composition[hossz][5]+=t_n[j][k];     composition[hossz][6]+=r_n[j][k];     composition[hossz][7]+=q_n[j][k];     len=0; max=0;     for(l=1; l<=hossz; l++)     {      dat=sq_seq[j][k][l];      if(dat==4)      {       len++;      }      if(dat!=4)      {       if(len>max) max=len;       len=0;      }     }     if(hossz==1) max=len;     composition[hossz][8]+=max;    }    else    {     arany[0]+=1;    }   }  }  for(j=0; j<=MAXL_N; j++)  {   for(k=1; k<=8; k++)   {    if (arany[j]) composition[j][k]=composition[j][k]/arany[j];   }  }  for(j=0; j<=MAXL_N; j++)  {   arany[j]=arany[j]/SIZE;  }  tf=fopen("result2.dat","a+t");   for(j=1; j<=(MAXL_N-1); j++)  {   fprintf(tf,"%8.4f", composition[j][1]);  }  fprintf(tf,"%8.4f\n", composition[MAXL_N][1]);  for(j=1; j<=(MAXL_N-1); j++)  {   fprintf(tf,"%8.4f", composition[j][2]);  }  fprintf(tf,"%8.4f\n", composition[MAXL_N][2]);  for(j=1; j<=(MAXL_N-1); j++)  {   fprintf(tf,"%8.4f", composition[j][3]);  }  fprintf(tf,"%8.4f\n", composition[MAXL_N][3]);  for(j=1; j<=(MAXL_N-1); j++)  {   fprintf(tf,"%8.4f", composition[j][4]);  }  fprintf(tf,"%8.4f\n", composition[MAXL_N][4]);  for(j=1; j<=(MAXL_N-1); j++)  {   fprintf(tf,"%8.4f", composition[j][5]);  }  fprintf(tf,"%8.4f\n", composition[MAXL_N][5]);  for(j=1; j<=(MAXL_N-1); j++)  {   fprintf(tf,"%8.4f", composition[j][6]);  }  fprintf(tf,"%8.4f\n", composition[MAXL_N][6]);  for(j=1; j<=(MAXL_N-1); j++)  {   fprintf(tf,"%8.4f", composition[j][7]);  }  fprintf(tf,"%8.4f\n", composition[MAXL_N][7]);  for(j=1; j<=(MAXL_N-1); j++)  {   fprintf(tf,"%8.4f", composition[j][8]);  }  fprintf(tf,"%8.4f\n", composition[MAXL_N][8]);  fprintf(tf,"%8d", lepesszam);   for(j=0; j<=(MAXL_N-1); j++)  {   fprintf(tf,"%9.5f", arany[j]);  }  fprintf(tf,"%9.5f\n", arany[MAXL_N]);  fflush(tf);  fclose(tf); } void result3(void) {  tf=fopen("result3.dat","a+t");   for(j=1; j<=SIDE; j++)  {   for(k=1; k<=SIDE; k++)   {    for(l=1; l<=(sq_length[j][k]); l++)    {     fprintf(tf," %1d", sq_seq[j][k][l]);    }    fprintf(tf," 0\n");   }  }  fflush(tf);  fclose(tf); }  main()
 {
  inic();                               /* initialisation */
  for(i=0; i<=STEPNUMBER; i++)          /* Monte-Carlo cycle */
  {
   if((i%SAVEL)==0) {save();}   if((i%RESULTL1)==0) {result1();}   if((i%RESULTL2)==0) {result2();}   if((i%RESULTL3)==0) {result3();}   mc_step();
  }
 }