/* fourier.c  94.11.20
 * Copyright 1983-1992   Albert Davis
 * all functions needed in addition to the transient analysis
 * to perform the fourier command.
 */
#include "ecah.h"
#include "error.h"
#include "io.h"
#include "mode.h"
#include "probh.h"
#include "declare.h"
/*--------------------------------------------------------------------------*/
	void	foinit(long);
	void	fonext(void);
	void	fosave(double);
	void	foout(double,double,double);
static  void	fo_head(int);
static	polar_t	fo_convert(complex_t*,polar_t*,int,int);
static	void	fo_print(polar_t*,int,int,double,polar_t);
static  polar_t	topolar(complex_t);
static  int	stepno(double,double,double);
/*--------------------------------------------------------------------------*/
extern const struct ioctrl io;
static complex_t *fodata;	/* fourier data ..			    */
				/* set up as an array [timesteps][probes]   */
extern probe_t *probelist[];
extern char e_om[], e_int[];
static int probes;	/* number of probes (from probe cmd) */
static long timesteps;	/* number of time steps in tran analysis, incl 0 */
static long datapts;	/* total number of data points (timepts * probes) */
static long indx;
static int timestepno;
/*--------------------------------------------------------------------------*/
/* foinit:  initialize fourier analysis
 * allocate memory, set up arrays, etc.
 */
void foinit(long points)
{
 long statsize;
 
 qfree((void**)&fodata);
 timesteps = points + 1;
 probes = count_probes(probelist[sFOURIER]);
 datapts = probes * timesteps;
 statsize = datapts * sizeof(complex_t);
 if (statsize > SEGSIZ)
    error(bERROR, "too many data points\n");
 fodata = (complex_t*)calloc((unsigned int)statsize, 1);
 if (!fodata)
    error(bERROR,  e_om, "fourier");
 timestepno = 0;
 indx = 0;
}
/*--------------------------------------------------------------------------*/
void fonext(void)                         /* next time step		    */
{                                         /* set to next column in 2d array.*/
 indx = ++timestepno;	  		  /* stored by row		    */
}                                         /* ignore the first (time)	    */
/*--------------------------------------------------------------------------*/
/* fosave:  save single point of transient analysis
 * no need to save imag value.  always 0
 * all points at a probe are a vector
 */
void fosave(double realvalue)
{
 if (fodata){
    if (indx >= datapts)
	error(bERROR, e_int, "fourier");
    fodata[indx].x = realvalue;
    indx += timesteps;
 }
}
/*--------------------------------------------------------------------------*/
/* foout:  print out the results of the transform
 */
void foout(double sstart, double sstop, double sstep)
{
 int prob;			/* probe counter.  step thru list */
 int startstep, stopstep;	/* step number (0 and 9 for Spice defaults) */

 plclose();
 plclear();
 startstep = stepno(0., sstep, sstart);
 stopstep  = stepno(0., sstep, sstop );

 for (prob = 0;  prob < probes;  prob++){
    polar_t maxvalue;
    polar_t *pdata;
    complex_t *cdata;
    
    cdata = &fodata[prob*timesteps];
    pdata = (polar_t*)cdata;
    fo_head(prob);
    fft(cdata,  (int)timesteps-1,  0);
    maxvalue = fo_convert(cdata,pdata,startstep,stopstep);
    fo_print(pdata,startstep,stopstep,sstep,maxvalue);
 }
}
/*--------------------------------------------------------------------------*/
/* fo_head: print output header
 * arg is index into probe array, to select probe name
 */
static void fo_head(int prob)
{
 mprintf(io.where,"# %-10s", probename(probelist[sFOURIER][prob]));
 mprintf(io.where,"--------- actual ---------  -------- relative --------\n");
 mprintf(io.where,"#freq       ");
 mprintf(io.where,"value        dB      phase  value        dB      phase\n");
}
/*--------------------------------------------------------------------------*/
/* fo_convert: convert data array from complex to polar, (copy)
 * return max magnitude, for scaling
 */
static polar_t fo_convert(
	complex_t *cdata,
	polar_t *pdata,
	int startstep,
	int stopstep)
{
 polar_t maxvalue;
 int fstep;
 maxvalue.mag = maxvalue.phz = 0.;
 for (fstep = startstep;  fstep <= stopstep;  fstep++){
    pdata[fstep] = topolar(cdata[fstep]);
    if (fstep > 0)
       pdata[fstep].mag *= 2;
    if (pdata[fstep].mag > maxvalue.mag){
       maxvalue = pdata[fstep];
    }
 }
 if (maxvalue.mag == 0.){
    maxvalue.mag = 1.;
    maxvalue.phz = 0.;
 }
 return maxvalue;
}
/*--------------------------------------------------------------------------*/
/* fo_print: print results of fourier analysis
 * for all points at single probe
 */
static void fo_print(
	polar_t *data,		/* data array, polar coordinates */
	int startstep,		/* first step to print, usually 0 */
	int stopstep,		/* last step to print, # of harmonics */
	double sstep,		/* frequency step size */
	polar_t maxvalue)	/* scale factor for relative */
{
 int fstep;
 for (fstep = startstep;  fstep <= stopstep;  ++fstep){
    polar_t unscaled;
    polar_t scaled;
    double freq = sstep * fstep;

    unscaled = data[fstep];
    scaled.mag = unscaled.mag / maxvalue.mag;
    scaled.phz = unscaled.phz - maxvalue.phz;
    if (scaled.phz > 180.)
       scaled.phz -= 360.;
    if (scaled.phz < -180.)
       scaled.phz += 360.;
    mprintf(io.where, "%s%s%7.2f %8.3f %s%7.2f %8.3f\n",
       ftos(freq,        "           ",5,io.formaat),
       ftos(unscaled.mag,"           ",5,io.formaat),
       db(unscaled.mag),
       unscaled.phz,
       ftos(scaled.mag,  "           ",5,io.formaat),
       db(scaled.mag),
       scaled.phz ) ;
 }
}
/*--------------------------------------------------------------------------*/
/* topolar: convert rectangular to polar
 * rotates 90 degrees!  (ref to sine instead of cosine)
 */
static polar_t topolar(complex_t arg)
{
 polar_t rv;
 rv.mag = hypot( arg.x, arg.y);
 rv.phz = phase(-arg.y, arg.x);
 return rv;
}
/*--------------------------------------------------------------------------*/
/* stepno: return step number given its frequency or time
 */
static int stepno(double start, double step, double here)
{
 double num;
 num = (here-start)/step + .5;
 return (int)num;
}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
