/* SOMakeWaveBasis -- Builds a test basis of SOWaveFunctios. */
/* Last edited on 2003-07-04 17:22:38 by ra014852 */

#include <SOFunction.h>
#include <SOWaveFunction.h>
#include <SOParams.h>
#include <SOBasic.h>
#include <SOPlot.h>
#include <SOPlotParams.h>

#include <js.h>
#include <ps.h>
#include <stdio.h>
#include <math.h>
#include <values.h>

typedef struct Options
  { 
    char *outName;         /* Name of basis output file (minus the ".bas" extension) */
    double maxFreq;        /* Value used to control the range of frequences relative to {N} */    
    /* Plotting options: */
    PlotOptions plt;      /* Plot format and style options. */
  } Options;
  
/* INTERNAL PROTOTYPES */

Options *GetOptions(int argn, char **argc);  
SOPlot_PageState *SetPlot(PlotOptions *plt);
FILE *OpenPlot(char *name, SOPlot_PageState *pgs);

/* IMPLEMENTATIONS */

int main(int argn, char **argc)
  { int i, j, nel = 0;
    Options *o = GetOptions(argn, argc);
    int maxFI = (int)ceil(o->maxFreq);
    Basis bas = SOFunctionRef_vec_new(100);
    int k[2]; /* Frequency vector. */

    /* 2D Plotting */
    SOPlot_PageState *pgs = SetPlot(&(o->plt));
    FILE *psf = OpenPlot(txtcat(o->outName, "-freqs"), pgs);

    for(i = -maxFI; i <= +maxFI; i++)
      { for(j = -maxFI; j <= +maxFI; j++)
        { double f2 = i * i + j * j;
          k[0] = i; k[1] = j;
          if(sqrt(f2) <= o->maxFreq)
            { SOWaveFunction *wf = SOWaveFunction_FromFreq(2, k);
              nel++; 
	      SOFunctionRef_vec_expand(&bas, nel);
              bas.el[nel-1] = (SOFunction *)wf;

              /* 2D Plotting */
              { double xp = (double)i/((double)maxFI);
                double yp = (double)j/((double)maxFI);
                ps_fill_dot(psf, xp, yp,  0.5, 1, 0, 0);
              }
            }
        }
      } 
    SOFunctionRef_vec_trim(&bas, nel);

    FILE *wr = open_write(txtcat(o->outName, ".bas"));

    SOFunction_WriteBasis(wr, bas);

    fclose(wr);

    /* 2D Plotting */
    SOPlot_Close(psf, pgs);

    return 0;
  }
  

SOPlot_PageState *SetPlot(PlotOptions *plt)
  {
    double hFigSize = plt->figSize;             /* In mm. */
    double vFigSize = plt->figSize * SQRTHALF;  /* In mm. */ 
    double figMargin = (plt->eps ? 0.5 : 2.0);  /* In mm. */
    int captionLines = (plt->eps ? 0 : 1);      /* No caption. */
    int hCount = 1, vCount = 1;  /* Figures per row and column. */
    SOPlot_PageState *pgs;

    pgs = SOPlot_PageState_New
      ( plt->eps, 
        plt->paperSize, hFigSize, vFigSize, figMargin, figMargin,
        captionLines, hCount, vCount
      );

    return pgs;
  }


FILE *OpenPlot(char *name, SOPlot_PageState *pgs)
  {
    double minR[2], maxR[2];     /* Rectangle to plot. */
    SOPlot_File *psf = NULL;

    minR[X] = -1.05; maxR[X] = +1.05;
    minR[Y] = -1.05; maxR[Y] = +1.05;
   
    psf = SOPlot_BeginFigure(
      psf, pgs, name, name, 
      minR[X], maxR[X], 
      minR[Y], maxR[Y]
    );
  
    ps_draw_segment(psf, minR[X], 0.0, maxR[X], 0.0);
    ps_draw_segment(psf, 0.0, minR[Y], 0.0, maxR[Y]);
    ps_draw_rectangle(psf, minR[X], maxR[X], minR[Y], maxR[Y]);

    return psf;
  }
  

#define PPUSAGE SOParams_SetUsage

Options *GetOptions(int argn, char **argc)
  {
    Options *o = (Options *)notnull(malloc(sizeof(Options)), "no mem");
    SOParams_T *pp = SOParams_NewT(stderr, argn, argc);
    PPUSAGE(pp, "SOMakeWaveBasis \\\n");
    PPUSAGE(pp, "  -outName BASISNAME\n");
    PPUSAGE(pp, "  -maxFreq CONTROL_VAL\n");
    PPUSAGE(pp, SOPlotParams_FunctionHelp " \n");

    SOParams_GetKeyword(pp, "-outName");                               
    o->outName = SOParams_GetNext(pp);  

    SOParams_GetKeyword(pp, "-maxFreq");                               
    o->maxFreq = SOParams_GetNextDouble(pp, 0, MAXDOUBLE); 

    /* Plotting options: */
    o->plt = SOPlotParams_FunctionParse(pp);

    SOParams_Finish(pp);
    return o;
  }