/* @(#)fibrand.c	17.1.1.1 (ESO-DMD) 01/25/02 17:11:16 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Correspondence concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT    (c)  1988   European Southern Observatory
.IDENT        fibrand.c
.LANGUAGE     C
.AUTHOR       K. Banse   ESO/IPG
.KEYWORDS     Bulk data frames, Images
.PURPOSE      rearrange planes of a cube
.VERSION      1.00  950921:   Creation 
---------------------------------------------------------------------*/

#define NINT(a)      ((a) < 0 ? (int)((a) - 0.5 ) : (int)((a) + 0.5))

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



/*

*/

/***********************************************************************
*
*   Lagged Fibonnacci random number generator 
*
*   Paul Coddington, April 1993.
*
*
* General lagged Fibonnacci generator using subtraction, with lags 
* p and q, i.e. F(p,q,-) in Marsaglia's notation.
* 
* The random numbers X_{i} are obtained from the sequence:
* 
*    X_{i} = X_{i-q} - X_{i-p}   mod M
*  
* where M is 1 if the X's are taken to be floating point reals in [0,1),
* as they are here.
* 
* For good results, the biggest lag should be at least 1000, and probably
* on the order of 10000.
* 
* The following lags give the maximal period of the generator, which is 
* (2^{p} - 1)2^{n-1} on integers mod 2^n or reals with n bits in the 
* mantissa (see Knuth, or W. Zerler, Information and Control, 15 (1969) 67, 
* for a complete more list).
* 
*     P     Q
*   9689   471 
*   4423  1393
*   2281   715
*   1279   418
*    607   273
*    521   168
*    127    63
* 
* This program is based on the implementation of RANMAR in 
* F. James, "A Review of Pseudo-random Number Generators", 
* Comput. Phys. Comm. 60, 329 (1990).
* 
* For more details, see:
* 
* D.E. Knuth, The Art of Computer Programming Vol. 2: 
* Seminumerical Methods, (Addison-Wesley, Reading, Mass., 1981).
* 
* P. L'Ecuyer, Random numbers for simulation, Comm. ACM 33:10, 85 (1990).
* 
* G.A. Marsaglia, A current view of random number generators,
* in Computational Science and Statistics: The Interface,
* ed. L. Balliard (Elsevier, Amsterdam, 1985).
* 
***********************************************************************/
 


#define MAXSEED 900000000

#define SIGBITS 24   /* Number of significant bits */

/* Lags */
#define P 1279
#define Q  418

/* Variables for lagged Fibonacci */
double u[P+1];  /* seed table */
int pt0, pt1;   /* pointers into the seed table */



/***************************************************************************
	Initialize the random number generator.
        Taken from RMARIN in James's review -- initializes
	every significant bit using a combination linear 
	congruential and small lag Fibonacci generator.
***************************************************************************/
void RandomSeed(seed)
{
      int ij, kl, i, j, k, l;
      int ii, jj, m;
      double t, s;
      
      if (seed < 0) seed = - seed;
      seed = seed % MAXSEED;

      ij = seed / 30082;
      kl = seed - 30082 * ij;
      i = ((ij/177)% 177) + 2;
      j =  ij%177 + 2;
      k = ((kl/169)% 178) + 1;
      l =  kl%169;
      
      for ( ii = 1 ; ii <= P; ii++ ){
        s = 0.0;
        t = 0.5;
        for ( jj = 1 ; jj <= SIGBITS; jj++ ){
          m = (((i*j)% 179)*k)% 179;
          i = j;
          j = k;
          k = m;
          l = (53*l+1) % 169;
          if ( ((l*m)%64) >= 32) 
            s = s + t;
          t = 0.5 * t;
        }
        u[ii] = s;
      }
      pt0 = P;
      pt1 = Q;

      return;
}
  
/***************************************************************************
	Return a random double in [0.0, 1.0)
***************************************************************************/
double RandomNumber()
{
      double uni;
      
      uni = u[pt0] - u[pt1];
      if (uni < 0.0) uni = uni + 1.0;
      u[pt0] = uni;
      pt0 = pt0 - 1;
      if (pt0 == 0) pt0 = P;
      pt1 = pt1 - 1;
      if (pt1 == 0) pt1 = P;
            
      return(uni);
}

/*

*/


main()

{
char   *pntra;
char   ident[74], outfile[84], reffile[84], cbuf[12];

int  in[4], nobytes, nval, unit, null, fmtin;
int  naxis, npix[3], imnob, imnow, size, mapsize, pixoff;
long int seed;
register int  nr;

float  *fpntr, rr[4];

double  start[3], step[3], dst[6];




(void) SCSPRO("fibrand");
(void) SCKGETC("IN_A",1,80,&nval,reffile);		/* get ref frame */
(void) SCKGETC("OUT_A",1,80,&nval,outfile);             /* get output frame */
(void) SCKRDI("SEED",1,1,&nval,in,&unit,&null);
seed = (long int) in[0];

npix[1] = npix[2] = 1;

if (reffile[0] == '+')
   {
   (void) SCKRDI("ID",1,4,&nval,in,&unit,&null);
   (void) SCKRDD("RD",1,6,&nval,dst,&unit,&null);
   naxis = in[0];
   for (nr=0; nr<naxis; nr++)
      {
      npix[nr] = in[nr+1];
      start[nr] = dst[nr];
      step[nr] = dst[nr+naxis];
      }
   }
else
   {
   (void) SCFOPN(reffile,D_R4_FORMAT,0,F_IMA_TYPE,&imnob);
   (void) SCDRDI(imnob,"NAXIS",1,1,&nval,&naxis,&unit,&null);
   (void) SCDRDI(imnob,"NPIX",1,naxis,&nval,npix,&unit,&null);
   (void) SCDRDD(imnob,"START",1,naxis,&nval,start,&unit,&null);
   (void) SCDRDD(imnob,"STEP",1,naxis,&nval,step,&unit,&null);
   (void) SCFCLO(imnob);
   }


/* create frame */

size = npix[0] * npix[1] * npix[2];		/*  get total size of frame */
(void) SCFCRE(outfile,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,size,&imnob);
(void) SCDWRI(imnob,"NAXIS",&naxis,1,1,&unit);
(void) SCDWRI(imnob,"NPIX",npix,1,naxis,&unit);
(void) SCDWRD(imnob,"START",start,1,naxis,&unit);
(void) SCDWRD(imnob,"STEP",step,1,naxis,&unit);
if (naxis == 1)
   (void) strcpy(ident,"data          x-pix");
else if (naxis == 2)
   (void) strcpy(ident,"data          x-pix           y-pix ");
else
   (void) strcpy(ident,"data          x-pix           y-pix           z-pix ");

nval = (naxis + 1) * 16;
(void) SCDWRC(imnob,"CUNIT",1,ident,1,nval,&unit);
(void) strcpy(ident,"Random image (Fibonacci method) ");
(void) SCDWRC(imnob,"IDENT",1,ident,1,72,&unit);
rr[0] = rr[2] = rr[3] = 0.0;  rr[1] = 1.0;
(void) SCDWRR(imnob,"LHCUTS",rr,1,4,&unit);

RandomSeed(seed);			/* init Random Number Generator */


/* allocate virtual memory for buffer with right data types */


(void) SCKRDI("MONITPAR",20,1,&nval,&mapsize,&unit,&null);
mapsize *= mapsize;				/* square image */
mapsize *= 5;					/* cube of 5 planes also o.k. */
if (mapsize < size)
   {
   nval = mapsize * sizeof(float);
   pntra = malloc((unsigned long int)nval);
   if (pntra == (char *) 0) SCETER(33,"Could not allocate memory...");
   pixoff = 1;

loop:
   fpntr = (float *) pntra;
   for (nr=0; nr<mapsize; nr++)
      {
      *fpntr++ = (float)RandomNumber();
      }
   (void) SCFPUT(imnob,pixoff,mapsize,pntra);
   pixoff += mapsize;
   size -= mapsize;
   if (size >0)
      {
      if (size < mapsize) mapsize = size;
      goto loop;
      }
   }
else
   {
   (void) SCFMAP(imnob,F_O_MODE,1,size,&nval,&pntra);

   fpntr = (float *) pntra;
   for (nr=0; nr<size; nr++)
      {
      *fpntr++ = (float)RandomNumber();
      }
   }

(void) SCSEPI();
}