/*
 *         SINMAX  v0.2
 *
 * Funcao: cria um arquivo FITS de dimensao naxes[0] x naxes[1]
 *         contendo Nos maximos periodicos.
 *
 * Uso: basta invocar o nome e informar o numero de nos
 *
 * Dependencias: CFITSIO 
 *
 * Autor: Fabricio Ferrari, abr/maio 1999
 *
 * Historia: v0.1->v0.11: * ususario deve informar numero de nos.
 * 	     v0.11->v0.2: * Os maximos sao aleatoriamente escolhidos
 *			    num intervalo razoavel (keyword PICO)
 *			  * Agora cada 'estrela' e' calculada 
 *			    separadamente a partir de uma gaussiana
 *			    centrada em algum lugar dentro da caixa 
 *	  		    de tamanho dimensao/nos de maximo PICO
 */


#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include "fitsio.h"

#define SHIFT 0.2                         /* fraction of center displacement */
#define PICO 10                           /* maximo da gaussiana */
#define nrand() ((float) rand()/RAND_MAX) /* randomico ate 1     */
#define sq(x) (x*x)                       /* melhor que pow(x,2) */

struct {
  int   centrox;
  int   centroy;
  float sizex;
  float sizey;
  float sigmax;
  float sigmay;
  float maximo;
} quad;


int main(void){

  fitsfile *fptr;
  int status=0, i, j, u, v;
  long fpixel=1, naxis=2, nelements, exposure=1, airmass=1;
  long naxes[]={600,600};
  float array[600][600];
  unsigned int nos=8;

  fits_create_file( &fptr, "!sinmax.fits", &status);
  fits_create_img(fptr, FLOAT_IMG, naxis, naxes, &status);
  fits_update_key(fptr, TLONG, "EXPTIME", &exposure, "Total Exposure Time", &status);
  fits_update_key(fptr, TLONG, "AIRMASS", &airmass, "Air Mass", &status);
  

  printf("\n\n+-------- SinMax --------+\n");
  printf("\nQual o n�mero de m�ximos: ");
  scanf("%ud", &nos);


  quad.sizex = naxes[0]/nos;
  quad.sizey = naxes[1]/nos;
  printf("sizes  x:%f  y:%f\n", quad.sizex, quad.sizey);
  
  srandom(PICO); /* thought in something better as seed ?!...*/

  for(v=0; v<nos; v++){
    for(u=0; u<nos; u++){

      /* calcultate the square center plus a random displacement  */
      quad.centrox = quad.sizex*(1/2 + u + SHIFT*nrand());
      quad.centroy = quad.sizey*(1/2 + v + SHIFT*nrand());
      quad.sigmax  = 0.1 * quad.sizex*(1+nrand()); /* Sigma is greater than 0.1 box size */
      quad.sigmay  = 0.1 * quad.sizey*(1+nrand()); /* and smaller than 0.2 box size      */
      quad.maximo  = PICO * nrand();
      
      printf("u:%i v:%i cx:%3i cy:%3i sx:%02.1f sy:%02.1f max:%02.1f rand:%01.2g\n",\
	      u,v,quad.centrox,quad.centroy,quad.sigmax,quad.sigmay,quad.maximo,nrand());

      for (j=0; j< naxes[1]; j++){
	for (i=0; i<naxes[0]; i++){
	  array[j][i] = array[j][i] +  quad.maximo * \
	                  exp(-( sq((i - quad.centrox)/quad.sigmax) + \
			         sq((j - quad.centroy)/quad.sigmay) ));
	  /*  printf("%f ",array[j][i]); */
	} /* i loop end */ 
      } /* j loop end */
    } /* u loop end */
  } /* v loop end */


  nelements = naxes[0] * naxes[1];
  fits_write_img(fptr, TFLOAT, fpixel, nelements, array[0], &status);
  fits_close_file(fptr, &status);
  fits_report_error(stderr, status);

  return(status);

}