function myRRgen(tmax,Pe,Pn,seed)

 % tmax : maximum time
 % Pe : Probability of ectopy ~ 1 per hr 
 % Pn : Probability of noise ~ 16 per hr 
 % seed : a 32-bit random variable that can be used to initialize the generator 

 t = 0;	           % sum of intervals since the beginning of the simulation, in seconds 
 ts = 0;	        % t, in sample intervals */
 tsp = 0;	        % previous value of ts */
   
% Initialise and runn RR interval generator, pass in length and
% probability of artefacts and nosie 
initialize(seed, tmax, Pe, Pn);

% generate RR interval */
while ((t += generate()) < tmax) {	% add RR interval to running time */
    % calculate and output a quantized RR interval */
    ts = ()(SPS*t + 0.5);
    printf("%5.3f\n", (ts - tsp)/(()SPS));
    tsp = ts;
end


end % main





function new_rr = modify_ectopic( rr2,  rr1)
 % shorten beat by 80% +/- 10% randomly chosen */
  new_rr;
  delta;
  rrint;

  delta = (()rand()/RAND_MAX);
  delta = 0.2*(delta-0.5);
  % length of last rr interval */
  rrint = rr2 - rr1;
  % shorten it by 0.7 +/-  a bit */
  rrint = (0.7 - delta)*rrint; 
  % work out new timing */
  new_rr = rr1 + rrint;

end

%---------------------------------------------------------------------------*/
% generate an extra beat (artefact)                                         */
%---------------------------------------------------------------------------*/

 make_noise( rr2,  rr1)
{ % shorten beat by 80% +/- 10% randomly chosen */
   new_rr;
   delta;
   rrint;

  delta = (()rand()/RAND_MAX);
  delta = ((delta*(SPS-(2.0/SPS))) + (1.0/SPS))/SPS; 

  % length of last rr interval */
  rrint = rr2 - rr1;
  % shorten it by 0.5 +/-  a bit */
  rrint = delta*rrint; 
  % work out new timing */
  new_rr = rr1 + rrint;

  return(new_rr); 
}

%--------------------------------------------------------------------------*/
% CHECK IF ECTOPICS SHOULD BE ADDED                                        */ 
%--------------------------------------------------------------------------*/

int check_for_ectopic( Pe)
{
  %  Pe=0.0003;*/ % probability of ectopic (~1 per hour)*/

  % no state or hr dependency */
  int rtn=NO; 
   rand_no;

  % check to see if we really want ectopic beats */
  if(Pe==0) return(NO);

  % if we do ... */

  % pick a random number between 0 and 1*/
  rand_no =  (()rand()/RAND_MAX);

  if(rand_no<Pe)    rtn = YES;
  else              rtn = NO;

  return(rtn);
}

%---------------------------------------------------------------------------*/
%      CHECK IF ARTEFACT SHOULD BE ADDED                                    */
%---------------------------------------------------------------------------*/

int check_for_noise(int state,  hr,  Pn)
{
   rand_no;
  int   rtn =NO; 
   Pa0 = 0.0004; % probability of artefact state 0 - ASLEEP*/
   Pa1 = 0.0048; % probability of artefact state 1 - AWAKE */

  Pa1 = Pn; 
  Pa0 = Pn/12;  % Probability of noise in sleep is 12 times less than when awake */

  % check to see if we really want to add noise beats */
  if(Pn==0) return(NO);

  % if we do ... */

  % pick a random number */
  rand_no =  (()rand()/RAND_MAX);

  % if asleep, flag an artefact randomly (scaling using heart rate) */
  if(state==ASLEEP){
    if((rand_no*60/hr)<Pa0)    rtn = YES;
    else                       rtn = NO;
  }

  % if awake, flag an artefact randomly with higher prob*/
  if(state!=ASLEEP){
    if((rand_no*60/hr)<Pa1)    rtn = YES;
    else                       rtn = NO;
  }

  return(rtn);
}

%---------------------------------------------------------------------------*/
%      ALLOCATE MEMORY FOR VECTOR                                           */
%---------------------------------------------------------------------------*/
 
 *vector( nl,  nh)
% allocate a  vector with subscript range v[nl..nh] */
{
         *v;
 
        v=( *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof()));
        if (!v) printf("allocation failure in vector");
        return v-nl+NR_END;
}

%---------------------------------------------------------------------------*/
%      FREE MEMORY FOR VECTOR                                               */
%---------------------------------------------------------------------------*/
 
void free_vector( *v,  nl,  nh)
% free a  vector allocated with vector() */
{
        free((FREE_ARG) (v+nl-NR_END));
}

%---------------------------------------------------------------------------*/
%      MEAN CALCULATOR                                                      */
%---------------------------------------------------------------------------*/
 
 mean( *x, int n)
% n-by-1 vector, calculate mean */
{
        int j;
         add;
 
        add = 0.0;
        for(j=1;j<=n;j++)  add += x[j];
 
        return (add/n);
}

%---------------------------------------------------------------------------*/
%      STANDARD DEVIATION CALCULATOR                                        */
%---------------------------------------------------------------------------*/
 
 std( *x, int n)
% n-by-1 vector, calculate standard deviation */
{
        int j;
         add,mean,diff,total;
 
        add = 0.0;
        for(j=1;j<=n;j++)  add += x[j];
        mean = add/n;
 
        total = 0.0;
        for(j=1;j<=n;j++)
        {
           diff = x[j] - mean;
           total += diff*diff;
        }
 
        return (sqrt(total/(n-1)));
}

%--------------------------------------------------------------------------*/
%    UNIFORM DEVIATES                                                      */
%--------------------------------------------------------------------------*/

 ran1( *idum)
{
	int j;
	 k;
	static  iy=0;
	static  iv[NTAB];
	 temp;

	if (*idum <= 0 || !iy) {
		if (-(*idum) < 1) *idum=1;
		else *idum = -(*idum);
		for (j=NTAB+7;j>=0;j--) {
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if (*idum < 0) *idum += IM;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ;
	*idum=IA*(*idum-k*IQ)-IR*k;
	if (*idum < 0) *idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j] = *idum;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}

%--------------------------------------------------------------------------*/
%    GAUSSIAN DEVIATES                                                     */
%--------------------------------------------------------------------------*/

 gasdev( *idum)
{
	 ran1( *idum);
	static int iset=0;
	static  gset;
	 fac,rsq,v1,v2;

	if  (iset == 0) {
		do {
			v1=2.0*ran1(idum)-1.0;
			v2=2.0*ran1(idum)-1.0;
			rsq=v1*v1+v2*v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return gset;
	}
}


%--------------------------------------------------------------------------*/
%    FFT                                                                   */
%--------------------------------------------------------------------------*/

void four1( data[], unsigned  nn, int isign)
{
	unsigned  n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	 tempr,tempi;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(data[j],data[i]);
			SWAP(data[j+1],data[i+1]);
		}
		m=n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=mmax << 1;
		theta=isign*(6.28318530717959/mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {
			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*data[j]-wi*data[j+1];
				tempi=wr*data[j+1]+wi*data[j];
				data[j]=data[i]-tempr;
				data[j+1]=data[i+1]-tempi;
				data[i] += tempr;
				data[i+1] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}

%--------------------------------------------------------------------------*/
%    BLOCK LENGTH                                                          */
%--------------------------------------------------------------------------*/

int blocklength()
{
    L;
   double alpha,beta,u;

   alpha = 5466.8;
   beta = - 2.2;

   u = (double)rand()/RAND_MAX;
   while(u < 1e-3)  u = (double)rand()/RAND_MAX;
   
   L = (int)pow( u/alpha, 1.0/beta);
 
   return L;
}

%--------------------------------------------------------------------------*/
%    TRANS LENGTH                                                          */
%--------------------------------------------------------------------------*/

int translength()
{
   int L;
    u;

   u = ()rand()/RAND_MAX;
   L = 5 + 25*u;

   return L;
}

%--------------------------------------------------------------------------*/
%    INTERP                                                                */
%--------------------------------------------------------------------------*/

void interp( *y,  *x, int n, int r)
{
   int i,j;
    a;

   for(i=1;i<=n-1;i++)
   {
      for(j=1;j<=r;j++) 
      {
         a = (j-1)*1.0/r;
         y[(i-1)*r+j] = (1.0-a)*x[i] + a*x[i+1];
      }
   }
   

}

%--------------------------------------------------------------------------*/
%    GENERATE TR BLOCK                                                     */
%--------------------------------------------------------------------------*/

void trblock( *tr,  flo,  fhi, 
 flostd,  fhistd,  lfhfratio,  
 rrmean,  rrtrend,  rrstd, int nbeats)
{
   int i,j,n;
    c1,c2,w1,w2,sig1,sig2,xstd,ratio;
    sfrr,dtrr;
    df,dw1,dw2,*w,*Hw,*Sw,*ph0,*ph,*SwC,*x,*rr;

   sfrr = 1.0;
   dtrr = 1.0;
   n = (int)pow(2.0, ceil(log(1.2*nbeats*rrmean*sfrr)/log(2.0))); 
   if(n < 256) n = 256;

   w = vector(1,n);
   Hw = vector(1,n);
   Sw = vector(1,n);
   ph0 = vector(1,n/2-1);
   ph = vector(1,n);
   SwC = vector(1,2*n);
   x = vector(1,n);
   rr = vector(1,n*10);

   w1 = 2.0*PI*flo;
   w2 = 2.0*PI*fhi;

   c1 = 2.0*PI*flostd;
   c2 = 2.0*PI*fhistd;
   sig2 = 1.0;
   sig1 = lfhfratio;


   df = sfrr/n;
   for(i=1;i<=n;i++) w[i] = (i-1)*2.0*PI*df;
   for(i=1;i<=n;i++) 
   {
      dw1 = w[i]-w1;
      dw2 = w[i]-w2;
      Hw[i] = sig1*exp(-dw1*dw1/(2.0*c1*c1))/sqrt(2*PI*c1*c1) 
            + sig2*exp(-dw2*dw2/(2.0*c2*c2))/sqrt(2*PI*c2*c2); 
   }

   for(i=1;i<=n/2;i++) Sw[i] = (sfrr/2.0)*Hw[i];
   for(i=n/2+1;i<=n;i++) Sw[i] = (sfrr/2.0)*Hw[n-i+1];

   % randomise the phases */
   for(i=1;i<=n/2-1;i++) ph0[i] = 2.0*PI*()rand()/RAND_MAX;
   ph[1] = 0.0;
   for(i=1;i<=n/2-1;i++) ph[i+1] = ph0[i];
   ph[n/2+1] = 0.0;
   for(i=1;i<=n/2-1;i++) ph[n-i+1] = - ph0[i]; 



   % make complex spectrum */
   for(i=1;i<=n;i++) SwC[2*i-1] = Sw[i]*cos(ph[i]);
   for(i=1;i<=n;i++) SwC[2*i] = Sw[i]*sin(ph[i]);


   % calculate inverse fft */
   four1(SwC,n,-1);

   % extract real part */
   for(i=1;i<=n;i++) x[i] = (1.0/n)*SwC[2*i-1];


   xstd = std(x,n);
   ratio = rrstd/xstd; 

   for(i=1;i<=n;i++) x[i] *= ratio;
   for(i=1;i<=n;i++) x[i] += rrmean;

   % upsample */
   interp(rr, x, n, 10);

   % incorporate trend */
   for(i=1;i<=n*10;i++) rr[i] = rr[i] + 0.1*rrtrend*(-0.5 + 1.0*i/(n*10));

   % add noise */
   for(i=1;i<=n*10;i++) rr[i] += (rr[i]/10.0)*()rand()/RAND_MAX;

   % determine r-peak times */
   tr[1] = rr[1];
   for(i=2;i<=nbeats;i++)
   {
      j = (int)ceil( (tr[i-1] + rr[i])/0.1); 
      tr[i] = tr[i-1] + rr[j];
   }


   free_vector(w,1,n);
   free_vector(Hw,1,n);
   free_vector(Sw,1,n);
   free_vector(ph0,1,n/2-1);
   free_vector(ph,1,n);
   free_vector(SwC,1,2*n);
   free_vector(x,1,n);
   free_vector(rr,1,n*10);
}


%--------------------------------------------------------------------------*/
%    GENERATE TRANSIENT                                                    */
%--------------------------------------------------------------------------*/

void trtrans( *tr, int n,  rr1,  rr2,  *trb)
{
   int i,iturn;
    rr12,rrthrs,rrdrop,rrmin,a,b,*rr,*rrb,*rrdev,rrbmean;

   rr = vector(1,n);
   rrb = vector(1,n);
   rrdev = vector(1,n);

   rr12 = MIN(rr1,rr2);


   rrthrs = 0.4 + 0.2*()rand()/RAND_MAX;
   rrdrop = 0.1 + 0.1*()rand()/RAND_MAX;
   rrmin = MAX(rrthrs, rr12-rrdrop); 

   % choose turning point at U(0.3,0.6) x (n+1) */
   iturn = (int)ceil((0.3 + 0.3*()rand()/RAND_MAX)*(n+1)); 
   
   % create triangle using slopes a and b */
   a = (rrmin - rr1)/(iturn - 0.0);
   b = (rr2 - rrmin)/(n+1 - iturn);

   for(i=1;i<=iturn;i++) rr[i] = rr1 + a*i;
   for(i=iturn+1;i<=n;i++) rr[i] = rrmin + b*(i-iturn);

   % make rr dynamics from block mean zero */
   rrb[1] = trb[1];
   for(i=2;i<=n;i++) rrb[i] = trb[i] - trb[i-1];
   rrbmean = mean(rrb, n);
   for(i=1;i<=n;i++) rrdev[i] = rrb[i] - rrbmean;

   tr[1] = rr[1] + rrdev[1];
   for(i=2;i<=n;i++) tr[i] = tr[i-1] + rr[i] + rrdev[i];

   free_vector(rr,1,n);
   free_vector(rrb,1,n);
   free_vector(rrdev,1,n);
}


% This function is called once, before generate() is called.  See main()
   for a description of its arguments.
*/

void initialize( seed,  tmax,  Pe,  Pn)
{
  %  srand((unsigned int)seed); 
      return; */
   
    i,j,time,Nb,Nt,Ndays,Nbeats;
    RRmean,RRamp,dRRamp,RRstdmin,RRstdmax,RRtrendamp,LFHFmin,LFHFmax;
    Tc,dRR,rrmean,rrstd,rrtrend,lfhfratio,rra,rrb;
    Tsleepstart,Tsleepdur,Tsleepfinish;
    flo,fhi,flostd,fhistd,sfrr,dtrr;
    *RR0,*tr0,*tr1,*trb,*trpeaks,u,RRsleep,acir;
    T1, T2, T3, T4;

    hr, *state;
   int   test;

   %   printf("seed = %ld\n",seed); */
   %   printf("tmax = %ld\n",tmax); */

   % set size in seconds and intialise global vector of r peaks */
   N = tmax;
   trpeaks_new = vector(1,2*N);


   % this memory allocation works for total averaged HR <= 120 bpm */
   RR0 = vector(1,2*N);
   trpeaks = vector(1,2*N);
   tr0 = vector(1,2*N);
   tr1 = vector(1,2*N);
   trb = vector(1,2*N);
   state = vector(1,2*N);

   rseed = -seed;

   % define RR  term characteristics */
   RRmean = 0.7 + 0.3*()rand()/RAND_MAX;
   RRamp = 0.075 + 0.2*()rand()/RAND_MAX;
   RRsleep = 0.1 + 0.1*()rand()/RAND_MAX;
   dRRamp = 0.03 + 0.1*()rand()/RAND_MAX;
   RRstdmin = 0.01;
   RRstdmax = 0.02;
   RRtrendamp = 1.0 + 0.25*()rand()/RAND_MAX;
   LFHFmin = 0.5;
   LFHFmax = 8.0;

    % define characteristics of rr power spectrum */
   flo = 0.1;
   fhi = 0.25;
   flostd = 0.01;
   fhistd = 0.01;
   sfrr = 1;
   dtrr = 1/sfrr;


   % define mean RR-interval using circadium rhythm period */
   Tc = (24.0 + gasdev(&rseed))*3600.0;
   acir = ()rand()/RAND_MAX;
   T1 = 10.0 + 3.0*()rand()/RAND_MAX;
   T2 = 7.0 + 2.0*()rand()/RAND_MAX;
   T3 = 4.0 + 2.0*()rand()/RAND_MAX;
   T4 = 2.0 + 2.0*()rand()/RAND_MAX;
   for(i=1;i<=N;i++){

     RR0[i] = RRmean + RRamp*sin(PI + (2.0*PI/Tc)*i) 
              + 0.25*acir*RRamp*sin((2.0*PI/(Tc/T1))*i) 
              + 0.25*RRamp*(1.0 - acir)*cos((2.0*PI/(Tc/T2))*i)
              + 0.25*acir*RRamp*sin((2.0*PI/(Tc/T3))*i) 
              + 0.25*RRamp*(1.0 - acir)*cos((2.0*PI/(Tc/T4))*i);

     state[i]=AWAKE;
   }


   % find number of 24 hour periods and insert sleep stages */
   Ndays = ()ceil(N/(24*3600));
   % printf("Ndays = %ld\n",Ndays); */
   for(j=1;j<=Ndays;j++)
   {

      % define sleep state Tsleepstart ~ U(14,16) Tsleepdur ~ U(6,8) */
      Tsleepstart = (j-1)*24*3600 + ()((14.0 + 2.0*()rand()/RAND_MAX)*3600.0);
      Tsleepdur = ()((6.0 + 2.0*()rand()/RAND_MAX)*3600.0);    
      Tsleepfinish = Tsleepstart + Tsleepdur;
      Tsleepfinish = MIN(Tsleepfinish, N);
      
      %     printf("Tsleepstart = %ld\n",Tsleepstart);  */
      %      printf("Tsleepfinish = %ld\n",Tsleepfinish); */
      for(i=Tsleepstart;i<=Tsleepfinish;i++){
	RR0[i] = RRmean + 0.5*RRsleep*(1.0+sin(2*PI/(20*60) )) ;
	% set sleep state */
	state[i] = ASLEEP;
      }
   }


   % add noise to mean RR value */
   for(i=1;i<=N;i++) RR0[i] += 0.2*RRamp*gasdev(&rseed);


   % initialise count in time and beats */
   time = 1;
   % i indexes beats rather than seconds now ! */
   i = 1;
 
   % generate first block */
   Nb = blocklength();


   % choose rr generation states: rrmean, rrstd, lfhfratio */
   % deviation from RR0 is given by dRR ~ U(0, dRRamp); */
   % dRR = dRRamp*()(2.0*rand()/RAND_MAX - 1.0); */
   
   % mimic segment switching in PRL 87, p168105 */
   u = ()rand()/RAND_MAX; 
   dRR = 0.75*dRRamp*(1.0 + exp(gasdev(&rseed))/10.0)*u/fabs(u);
   while(fabs(dRR) > dRRamp)
   {
      u = ()rand()/RAND_MAX;
      dRR = 0.5*dRRamp*(1.0 + ()exp(gasdev(&rseed))/10.0)*u/fabs(u);
   }
      
   rrmean = RR0[time] + dRR;
   rrtrend = RRtrendamp*(2*()rand()/RAND_MAX - 1); 
   rrstd = RRstdmin + (RRstdmax - RRstdmin)*()rand()/RAND_MAX;
   lfhfratio = LFHFmin + (LFHFmax-LFHFmin)*()rand()/RAND_MAX;


   % generate RR over this block */
   trblock(tr0, flo, fhi, flostd, fhistd, lfhfratio, rrmean, rrtrend, rrstd, Nb);

   for(j=1;j<=Nb;j++) trpeaks[i+j] = trpeaks[i] + tr0[j];
   i += Nb; 

   % convert beat time stamp back to nearest second */
   time = ()(trpeaks[i]); 
   rra = tr0[Nb] - tr0[Nb-1];

   % printf("time = %ld\n",time); */

   %
   for(j=2;j<=Nb;j++) printf("%d %e\n",j,tr0[j]- tr0[j-1]);
   exit(1); 
   */

   % generate record */
   while(time < N)
   {
    
      fhi = 0.2 + 0.1*()rand()/RAND_MAX;
      
      % choose a new block time in beats */
      Nb = blocklength();
      
      % choose a transition duration in beats */
      Nt = translength();

      % choose rr generation states: rrmean, rrstd, lfhfratio */
      dRR = dRRamp*()rand()/RAND_MAX;
      rrmean = RR0[time] + dRR;
      rrtrend = RRtrendamp*(2.0*()rand()/RAND_MAX - 1.0);
      rrstd = RRstdmin + (RRstdmax - RRstdmin)*()rand()/RAND_MAX;
      lfhfratio = LFHFmin + (LFHFmax - LFHFmin)*()rand()/RAND_MAX;

      % generate RR over this block */
    trblock(tr0, flo, fhi, flostd, fhistd, lfhfratio, rrmean, rrtrend, rrstd, Nb);
      rrb = tr0[2] - tr0[1];

     % generate RR during transition between blcoks with rr-intervals rra, rrb */
      rrtrend = 0.0;
    trblock(trb, flo, fhi, flostd, fhistd, lfhfratio, rrmean, rrtrend, rrstd, Nt);
    trtrans(tr1, Nt, rra, rrb, trb);


     % fill in transition and new block */
     for(j=1;j<=Nt;j++) trpeaks[i+j] = trpeaks[i] + tr1[j];
     i += Nt;
     for(j=1;j<=Nb;j++) trpeaks[i+j] = trpeaks[i] + tr0[j];
     i += Nb;

     % remember last rr-interval */
     rra = rrb;

     time = ()(trpeaks[i]);

   }
   Nbeats = i;


   % check for adding ectopics or artefacts */
   % loop over entire data set */
   j=1;
   trpeaks_new[1] = trpeaks[1];
   for(i=2;i<=Nbeats;i++){
     % find out what heart rate is at a particular rpeak time */
     hr=60.0/RR0[()(trpeaks[i])];     

     % check so see if we (randomly) want to make this ectopic */
     test = check_for_ectopic(Pe); % no state dependency */
     if (test!=NO){
       trpeaks_new[j] = modify_ectopic(trpeaks[i],trpeaks[i-1]);
     }    
     if (test == NO){ % just keep the current beat */
       test = check_for_noise(state[i], hr, Pn); % hr and state dependent */
       if (test!=NO){ % we have noise, add an extra beat */
	 trpeaks_new[j] = make_noise(trpeaks[i],trpeaks[i-1]);
	 j++;
       }
       % regardless, keep the current beat */
       trpeaks_new[j] = trpeaks[i];
     }
     j++;
   }

   free_vector(RR0,1,2*N);
   free_vector(trpeaks,1,2*N);
   free_vector(tr0,1,2*N);
   free_vector(tr1,1,2*N);
   free_vector(trb,1,2*N);
}


% This function is called once per RR interval.  It should return the length
   of the next (simulated) RR interval in seconds.
*/

 generate(void)
{
     rr;
    
    count++;
    rr= trpeaks_new[count+1]-trpeaks_new[count];

    % trap rounding errors */
    if(rr<(2.0/SPS)) rr = 2.0/SPS;

    return(rr);
}
    
int main(int argc, char **argv)
{
     t = 0.0;	    % sum of intervals since the beginning of the
			       simulation, in seconds */
     ts = 0;	    % t, in sample intervals */
     tsp = 0;	    % previous value of ts */
     tmax = 24*60*60;   % 24 hours, in seconds */
     seed;		    % a 32-bit random variable that can be used to
			       initialize the generator */
     Pe   = 0.0003;    % Probability of ectopy ~ 1 per hr */
     Pn   = 0.0048;    % Probability of noise ~ 16 per hr */
     atol();

    if (argc < 2) {
	fprintf(stderr, "usage: %s seed [tmax] P(ectopy) P(noise)\n", argv[0]);
	exit(1);
    }
    seed = atol(argv[1]);

    if (argc > 2)
	tmax = atol(argv[2]);
    else 
      tmax=115;

    if(tmax<115) 
        tmax=115;

    if(argc>3) Pe = atof(argv[3]);

    if(argc>4) Pn = atof(argv[4]);

    % Initialise and runn RR interval generator, pass in length and
       probability of artefacts and nosie */
    initialize(seed, tmax, Pe, Pn);

    % generate RR interval */
    while ((t += generate()) < tmax) {	% add RR interval to running time */
	% calculate and output a quantized RR interval */
	ts = ()(SPS*t + 0.5);
	printf("%5.3f\n", (ts - tsp)/(()SPS));
	tsp = ts;
    }

    exit(0);
}