#include<stdio.h>
#include<stdlib.h>
#include<time.h>

#define F77NAME(x) x## _
#define F77_NAME(name) name## __

struct linearad_common {
  double statelft,statergt,alpha,beta,ca,la,cr,lr,lc;
};
extern struct linearad_common F77NAME(linearad);

/*******************************************************************
 * This subroutine sets the initial conditions for the simulation. *
 *******************************************************************/
void reinitiate(int ncells, int ifirst, int ilast, double *xold)
{
  int col=(int)sqrt((double)ncells);
  double left=F77NAME(linearad).statelft;
  double right=F77NAME(linearad).statergt;
  double interv=(right-left)/(double)col;
  double speed, angle;
  int ic, count=0;
  double y_pos, x_pos;
  time_t unix_time;

  /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    Define the initial speed=sqrt(alpha/beta)
    But if beta=0, initial speed is set to 0
    +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
  if(F77NAME(linearad).beta!=0.)
    speed = sqrt(F77NAME(linearad).alpha/F77NAME(linearad).beta);
  else
    speed = 0.0;

  /*++++++++++++++++++++++++++++++++++++++++++++++++++
    Set initial positions and velocities:
    (Also initiate the random number generator)
  ++++++++++++++++++++++++++++++++++++++++++++++++++++*/
  y_pos = left;
  x_pos = left;
  srand(time(&unix_time));
  printf("time is %ld\n", unix_time);
  for(ic=ifirst; ic<=ilast; ic++) {
    /*+++++++++++++++++++++++++++++++++++
      Initial positions are initiated on
      a lattice square.
    ++++++++++++++++++++++++++++++++++++*/
    if(count>col) {
      count=0;
      y_pos+=interv;
      x_pos = left;
    }
    x_pos += interv;
    xold[ic] = -x_pos;
    xold[ic+ncells] = y_pos;
    count++;
    /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      Or we can set initial positions on a randomly distributed disc:
    +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
    /*
    xold[ic] = 100.;
    xold[ic+ncells] = 100.;
    for(;sqrt(xold[ic]*xold[ic]+xold[ic+ncells]*xold[ic+ncells])>1.0;) {
      xold[ic] = (double)rand()/(double)RAND_MAX*2.-1.;
      xold[ic+ncells] = (double)rand()/(double)RAND_MAX*2.-1.;
    }
    */


    /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      Initial velocities are initiated in random directions with
      the speed equal to sqrt(alpha/beta):
    ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
    /*
    angle = (double)rand()/(double)RAND_MAX*2*3.1415926535;
    xold[ic+2*ncells] = -speed*cos(angle);
    xold[ic+3*ncells] = speed*sin(angle);;
    */
    /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      Or we can set the initial velocities to zero:
    +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
    /*
    xold[ic+2*ncells] = 0.;
    xold[ic+3*ncells] = 0.;
    */
    /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      Or we can set the initial velocities in rotation:
    ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
    xold[ic+2*ncells] = - speed * xold[ic+ncells]
      / sqrt(xold[ic]*xold[ic]+xold[ic+ncells]*xold[ic+ncells]);
    xold[ic+3*ncells] = speed * xold[ic]
      / sqrt(xold[ic]*xold[ic]+xold[ic+ncells]*xold[ic+ncells]);
    /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      Or we can set the initial velocities to in radiation:
    ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
    /*
    xold[ic+2*ncells] = speed * xold[ic]
      / sqrt(xold[ic]*xold[ic]+xold[ic+ncells]*xold[ic+ncells]);
    xold[ic+3*ncells] = speed * xold[ic+ncells]
      / sqrt(xold[ic]*xold[ic]+xold[ic+ncells]*xold[ic+ncells]);
    */
  }
}