/* Steven Andrews, 9/9/98	*/
/* See documentation called Spectra doc */
/* Copyright 2003 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "Spectra.h"
#include "math2.h"
#include "Rn.h"
#include "Cn.h"
#include "Plot.h"
#include "DiskIO.h"
#include "RnSort.h"
#include "random.h"

sptr SpectAlloc(char *name,char *file,char *desc,char *xunit,char *yunit) {
	sptr s;
	
	s=(sptr) malloc(sizeof(struct spect));
	if(!s) return NULL;
	if(name)	strncpy(s->name,name,STRCHAR);
		else s->name[0]='\0';
	if(file)	strncpy(s->file,file,STRCHAR);
		else s->file[0]='\0';
	if(desc)	strncpy(s->desc,desc,STRCHAR);
		else s->desc[0]='\0';
	if(xunit)	strncpy(s->xunit,xunit,STRCHAR);
		else s->xunit[0]='\0';
	if(yunit)	strncpy(s->yunit,yunit,STRCHAR);
		else s->yunit[0]='\0';
	strncpy(s->color,"02",STRCHAR);
	s->n=0;
	s->cmplx=0;
	s->x=NULL;
	s->y=NULL;
	return s; }

void SpectFree(sptr s) {
	if(!s) return;
	if(s->x) free(s->x);
	if(s->y) free(s->y);
	free(s);
	return; }

int LoadSpect(sptr *spt,char *fname,int xcol,int ycol,int skip) {
	sptr s;
	float *a;
	int n,er;
	
	if(*spt) SpectFree(*spt);
	*spt=NULL;
	if(xcol<1 || ycol<1 || skip<0) return 6;
	s=SpectAlloc(fname,fname,fname,"x unit","y unit");
	if(!s) return 6;
	er=LoadData3(&a,&n,xcol,ycol,fname,skip);
	if(er) {SpectFree(s);return er;}
	s->x=allocV(n);
	s->y=allocV(n);
	if(!s->x || !s->y) {freeM(a);SpectFree(s);return 1;}
	s->n=n;
	columnM(a,s->x,n,2,0);
	columnM(a,s->y,n,2,1);
	sortV(s->x,s->y,n);
	freeM(a);
	*spt=s;
	return 0; }

int SaveSpect(sptr s) {
	float *a;
	int i,col,er;

	if(s->cmplx) col=3;
	else col=2;
	a=allocM(s->n,col);
	if(!a) return 1;
	for(i=0;i<s->n;i++) a[col*i]=s->x[i];
	if(s->cmplx)
		for(i=0;i<s->n;i++) {a[col*i+1]=s->y[2*i];a[col*i+2]=s->y[2*i+1];}
	else for(i=0;i<s->n;i++) a[col*i+1]=s->y[i];
	er=SaveData(a,s->n,col,s->file,0);
	if(er) s->file[0]='\0';
	freeM(a);
	return er; }

void SpectRange(sptr s,float *xa,float *xb,float *ya,float *yb,int fn) {
	float z;
	
	z=minV(s->x,s->n);
	if(!fn || (fn<0 && *xa<z) || (fn>0 && *xa>z)) *xa=z;
	z=maxV(s->x,s->n);
	if(!fn || (fn<0 && *xb>z) || (fn>0 && *xb<z)) *xb=z;
	z=minV(s->y,s->cmplx?2*s->n:s->n);
	if(!fn || (fn<0 && *ya<z) || (fn>0 && *ya>z)) *ya=z;
	z=maxV(s->y,s->cmplx?2*s->n:s->n);
	if(!fn || (fn<0 && *yb>z) || (fn>0 && *yb<z)) *yb=z;
	return; }

void TypeSpect(sptr s) {
	int n;
	float *x,*y;
	
	if(!s) {
		printf("NULL spectrum\n");
		return; }
	n=s->n;
	x=s->x;
	y=s->y;
	printf("name:%s\tfile:%s\tcolor:%s\n",s->name,s->file,s->color);
	printf("points:%i\txunit:%s\tyunit:%s\n",n,s->xunit,s->yunit);
	printf("desc:%s\n",s->desc);
	if(!x || !y)
		printf("No data.\n");
	else if(s->cmplx) {
		printf("real: (%f,%f)...(%f,%f)\n",x[0],y[0],x[n-1],y[2*n-2]);
		printf("imag: (%f,%f)...(%f,%f)\n",x[0],y[1],x[n-1],y[2*n-1]); }
	else printf("data: (%f,%f)...(%f,%f)\n",x[0],y[0],x[n-1],y[n-1]);
	return; }

void PlotSpect(sptr s) {
	int i;
	
	if(!s || !s->x || !s->y) return;
	SetColor(s->color[0]);
	PlotPt(s->x[0],s->y[0]);
	if(s->cmplx) {
		for(i=1;i<s->n;i++) PlotLine(s->x[i],s->y[2*i]);
		SetColor(s->color[1]);
		PlotPt(s->x[0],s->y[1]);
		for(i=1;i<s->n;i++) PlotLine(s->x[i],s->y[2*i+1]); }
	else for(i=1;i<s->n;i++) PlotLine(s->x[i],s->y[i]);
	SetColor('0'); }

int SpectMath(sptr s1,sptr s2,sptr *s3ptr,char *fn,float k) {
	sptr s3;
	int n,n2,cmplx,er,i,j,i2,sort;
	float *x1,*x2,*x3,*y1,*y2,*y3;
	float f1,f2,f3,*v1,*v2,num[MAXARG];
	sptr sz1,sz2;

	er=0;
	sort=0;
	if(!s1 || !s3ptr || !fn) return 6;
	if(*s3ptr) {
		s3=*s3ptr;
		if(s3->n!=s1->n || s3->cmplx!=s1->cmplx) {
			x3=allocV(s1->n);
			y3=allocV(s1->cmplx?2*s1->n:s1->n);
			if(!x3 || !y3) {er=1;goto failure;}
			if(s3->x) freeV(s3->x);
			if(s3->y) freeV(s3->y);
			s3->x=x3;
			s3->y=y3; }
		strncpy(s3->desc,s1->desc,STRCHAR);
		strncpy(s3->xunit,s1->xunit,STRCHAR);
		strncpy(s3->yunit,s1->yunit,STRCHAR); }
	else {
		s3=SpectAlloc(NULL,NULL,s1->desc,s1->xunit,s1->yunit);
		if(!s3) return 1;
		s3->x=allocV(s1->n);
		s3->y=allocV(s1->cmplx?2*s1->n:s1->n);
		if(!s3->x || !s3->y) {er=1;goto failure;}}
	n=s3->n=s1->n;
	cmplx=s3->cmplx=s1->cmplx;
	strncpy(s3->color,s1->color,STRCHAR);
	copyV(s1->x,s3->x,n);
	strncat(s3->desc," ",STRCHAR-strlen(s3->desc));
	strncat(s3->desc,fn,STRCHAR-strlen(s3->desc));
	x1=s1->x;
	x3=s3->x;
	y1=s1->y;
	y3=s3->y;
	if(s2) {
		x2=s2->x;
		y2=s2->y;
		n2=s2->n; }
	else {
		x2=y2=NULL;
		n2=0; }

	if(!strcmp(fn,"copy"))
		if(cmplx) copyV(y1,y3,2*n);
		else copyV(y1,y3,n);
	else if(!strcmp(fn,"smooth")) {
		if(cmplx) {er=8;goto failure;}
		else if(!isevenspV(x1,n,0.001)) {er=15;goto failure;}
		else if(k<0) {er=44;goto failure;}
		else if(!smoothV(y1,y3,n,(int)k)) {er=1;goto failure;}}
	else if(!strcmp(fn,"log"))
		if(cmplx) for(i=0;i<n;i++) {
				f1=CMPXmag(y1[2*i],y1[2*i+1]);
				y3[2*i]=(f1>0)?log10(f1):0;
				y3[2*i+1]=CMPXang(y1[2*i],y1[2*i+1])/2.302585093; }
		else for(i=0;i<n;i++)	y3[i]=(y1[i]>0)?log10(y1[i]):0;
	else if(!strcmp(fn,"exp"))
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=exp(y1[2*i])*cos(y1[2*i+1]);
				y3[2*i+1]=exp(y1[2*i])*sin(y1[2*i+1]); }
		else for(i=0;i<n;i++)	y3[i]=exp(y1[i]);
	else if(!strcmp(fn,"ln"))
		if(cmplx) for(i=0;i<n;i++) {
				f1=CMPXmag(y1[2*i],y1[2*i+1]);
				y3[2*i]=(f1>0)?log(f1):0;
				y3[2*i+1]=CMPXang(y1[2*i],y1[2*i+1]); }
		else for(i=0;i<n;i++)	y3[i]=(y1[i]>0)?log(y1[i]):0;
	else if(!strcmp(fn,"10^s")) {
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=exp(2.302585093*y1[2*i])*cos(2.302585093*y1[2*i+1]);
				y3[2*i+1]=exp(2.302585093*y1[2*i])*sin(2.302585093*y1[2*i+1]); }
		else for(i=0;i<n;i++)	y3[i]=pow(10,y1[i]); }
	else if(!strcmp(fn,"sqrt"))
		if(cmplx) for(i=0;i<n;i++) {
				f1=sqrt(CMPXmag(y1[2*i],y1[2*i+1]));
				f2=CMPXang(y1[2*i],y1[2*i+1])/2;
				y3[2*i]=f1*cos(f2);
				y3[2*i+1]=f1*sin(f2); }
		else for(i=0;i<n;i++)	y3[i]=(y1[i]>0)?sqrt(y1[i]):0;
	else if(!strcmp(fn,"x*s") || !strcmp(fn,"s*x"))
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=y1[2*i]*x3[i];
				y3[2*i+1]=y1[2*i+1]*x3[i]; }
		else for(i=0;i<n;i++)	y3[i]=y1[i]*x3[i];
	else if(!strcmp(fn,"s/x"))
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=x3[i]?y1[2*i]/x3[i]:0;
				y3[2*i+1]=x3[i]?y1[2*i+1]/x3[i]:0; }
		else for(i=0;i<n;i++)	y3[i]=x3[i]?y1[i]/x3[i]:0;
	else if(!strcmp(fn,"deriv1") || !strcmp(fn,"deriv"))
		if(cmplx) {er=8;goto failure;}
		else {
			v1=allocV(n);
			if(!v1) {er=1;goto failure;}
			deriv1V(y1,y3,n);
			deriv1V(x1,v1,n);
			divV(y3,v1,y3,n);
			freeV(v1); }
	else if(!strcmp(fn,"integ"))
		if(cmplx) {er=8;goto failure;}
		else {
			v1=allocV(n);
			if(!v1) {er=1;goto failure;}
			integV(y1,y3,n);
			deriv1V(x1,v1,n);
			multV(y3,v1,y3,n);
			j=-2;
			f1=interpolate1(x1,y3,n,&j,k);
			addKV(-f1,y3,y3,n);
			freeV(v1); }
	else if(!strcmp(fn,"deriv2"))
		if(cmplx) {er=8;goto failure;}
		else {
			v1=allocV(n);
			if(!v1) {er=1;goto failure;}
			deriv2V(y1,y3,n);
			deriv1V(x1,v1,n);
			multV(v1,v1,v1,n);
			divV(y3,v1,y3,n);
			freeV(v1); }
	else if(!strcmp(fn,"xderiv1"))
		if(cmplx) {er=8;goto failure;}
		else {
			v1=allocV(n);
			if(!v1) {er=1;goto failure;}
			divV(y1,x1,v1,n);
			deriv1V(v1,y3,n);
			deriv1V(x1,v1,n);
			divV(y3,v1,y3,n);
			multV(y3,x1,y3,n);
			freeV(v1); }
	else if(!strcmp(fn,"xderiv2"))
		if(cmplx) {er=8;goto failure;}
		else {
			v1=allocV(n);
			if(!v1) {er=1;goto failure;}
			divV(y1,x1,v1,n);
			deriv2V(v1,y3,n);
			deriv1V(x1,v1,n);
			multV(v1,v1,v1,n);
			divV(y3,v1,y3,n);
			multV(y3,x1,y3,n);
			freeV(v1); }
	else if(!strcmp(fn,"s+k") || !strcmp(fn,"k+s"))
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=y1[2*i]+k;
				y3[2*i+1]=y1[2*i+1]; }
		else for(i=0;i<n;i++)	y3[i]=y1[i]+k;
	else if(!strcmp(fn,"s-k"))
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=y1[2*i]-k;
				y3[2*i+1]=y1[2*i+1]; }
		else for(i=0;i<n;i++)	y3[i]=y1[i]-k;
	else if(!strcmp(fn,"k-s"))
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=k-y1[2*i];
				y3[2*i+1]=-y1[2*i+1]; }
		else for(i=0;i<n;i++)	y3[i]=k-y1[i];
	else if(!strcmp(fn,"s*k") || !strcmp(fn,"k*s"))
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=y1[2*i]*k;
				y3[2*i+1]=y1[2*i+1]*k; }
		else for(i=0;i<n;i++)	y3[i]=y1[i]*k;
	else if(!strcmp(fn,"k/s"))
		if(cmplx) for(i=0;i<n;i++) {
				f1=y1[2*i]*y1[2*i]+y1[2*i+1]*y1[2*i+1];
				y3[2*i]=f1?k*y1[2*i]/f1:0;
				y3[2*i+1]=f1?-k*y1[2*i+1]/f1:0; }
		else for(i=0;i<n;i++)	y3[i]=y1[i]?k/y1[i]:0;
	else if(!strcmp(fn,"s/k")) {
		if(!k) {er=9;goto failure;}
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=y1[2*i]/k;
				y3[2*i+1]=y1[2*i+1]/k; }
		else for(i=0;i<n;i++)	y3[i]=y1[i]/k; }
	else if(!strcmp(fn,"s^k"))
		if(cmplx) for(i=0;i<n;i++) {
				f1=CMPXmag(y1[2*i],y1[2*i+1]);
				f1=(f1==0 && k<=0)?0:pow(f1,k);
				f2=k*CMPXang(y1[2*i],y1[2*i+1]);
				y3[2*i]=f1*cos(f2);
				y3[2*i+1]=f1*sin(f2); }
		else for(i=0;i<n;i++)	y3[i]=((y1[i]==0 && k<=0) || (y1[i]<0 && !isinteger(k)))?0:pow(y1[i],k);
	else if(!strcmp(fn,"s*s"))
		if(!s2) {er=6;goto failure; }
		else if(cmplx && s2->cmplx)
			for(j=-2,i=0;i<n;i++) {
				f1=(x2[i]==x1[i])?y2[2*i]:interpolate1Cr(x2,y2,n2,&j,x1[i]);
				f2=(x2[i]==x1[i])?y2[2*i+1]:interpolate1Ci(x2,y2,n2,&j,x1[i]);
				y3[2*i]=y1[2*i]*f1-y1[2*i+1]*f2;
				y3[2*i+1]=y1[2*i]*f2+y1[2*i+1]*f1; }
		else if(!cmplx && !s2->cmplx)
			for(j=-2,i=0;i<n;i++)
				y3[i]=y1[i]*((x2[i]==x1[i])?y2[i]:interpolate1(x2,y2,n2,&j,x1[i]));
		else {er=10;goto failure;}
	else if(!strcmp(fn,"s/s"))
		if(!s2) {er=6;goto failure; }
		else if(cmplx && s2->cmplx)
			for(j=-2,i=0;i<n;i++) {
				f1=(x2[i]==x1[i])?y2[2*i]:interpolate1Cr(x2,y2,n2,&j,x1[i]);
				f2=(x2[i]==x1[i])?y2[2*i+1]:interpolate1Ci(x2,y2,n2,&j,x1[i]);
				f3=f1*f1+f2*f2;
				y3[2*i]=f3?(y1[2*i]*f1+y1[2*i+1]*f2)/f3:0;
				y3[2*i+1]=f3?(-y1[2*i]*f2+y1[2*i+1]*f1)/f3:0; }
		else if(!cmplx && !s2->cmplx)
			for(j=-2,i=0;i<n;i++) {
				f1=(x2[i]==x1[i])?y2[i]:interpolate1(x2,y2,n2,&j,x1[i]);
				y3[i]=f1?y1[i]/f1:0; }
		else {er=10;goto failure;}
	else if(!strcmp(fn,"s+s"))
		if(!s2) {er=6;goto failure; }
		else if(cmplx && s2->cmplx)
			for(j=-2,i=0;i<n;i++) {
				f1=(x2[i]==x1[i])?y2[2*i]:interpolate1Cr(x2,y2,n2,&j,x1[i]);
				f2=(x2[i]==x1[i])?y2[2*i+1]:interpolate1Ci(x2,y2,n2,&j,x1[i]);
				y3[2*i]=y1[2*i]+f1;
				y3[2*i+1]=y1[2*i]+f2; }
		else if(!cmplx && !s2->cmplx)
			for(j=-2,i=0;i<n;i++)
				y3[i]=y1[i]+((x2[i]==x1[i])?y2[i]:interpolate1(x2,y2,n2,&j,x1[i]));
		else {er=10;goto failure;}
	else if(!strcmp(fn,"s-s"))
		if(!s2) {er=6;goto failure; }
		else if(cmplx && s2->cmplx)
			for(j=-2,i=0;i<n;i++) {
				f1=(x2[i]==x1[i])?y2[2*i]:interpolate1Cr(x2,y2,n2,&j,x1[i]);
				f2=(x2[i]==x1[i])?y2[2*i+1]:interpolate1Ci(x2,y2,n2,&j,x1[i]);
				y3[2*i]=y1[2*i]-f1;
				y3[2*i+1]=y1[2*i]-f2; }
		else if(!cmplx && !s2->cmplx)
			for(j=-2,i=0;i<n;i++)
				y3[i]=y1[i]-((x2[i]==x1[i])?y2[i]:interpolate1(x2,y2,n2,&j,x1[i]));
		else {er=10;goto failure;}
	else if(!strcmp(fn,"merge"))
		if(!s2) {er=6;goto failure;}
		else if(cmplx && !s2->cmplx) {er=10;goto failure;}
		else if(!cmplx && s2->cmplx) {er=10;goto failure;}
		else {
			if(x2[0]<x1[0])	{
				sz1=s1;s1=s2;s2=sz1;
				x1=s1->x;x2=s2->x;
				y1=s1->y;y2=s2->y;
				n=s1->n;n2=s2->n; }
			i2=1+locateV(x2,x1[n-1],n2);
			if(i2==n2) {er=11;goto failure;}
			n+=n2-i2;
			x3=allocV(n);
			if(!x3) {er=1;goto failure;}
			y3=allocV(cmplx?2*n:n);
			if(!y3) {freeV(x3);er=1;goto failure;}
			freeV(s3->x);
			freeV(s3->y);
			s3->n=n;
			s3->x=x3;
			s3->y=y3;
			copyV(x1,x3,s1->n);
			copyV(x2+i2,x3+s1->n,n2-i2);
			if(cmplx) {
				for(i=0;x3[i]<x2[0];i++) {y3[2*i]=y1[2*i];y3[2*i+1]=y1[2*i+1];}
				for(j=-2;x3[i]<=x1[s1->n-1];i++) {
					y3[2*i]=y1[2*i]*(x2[i2]-x3[i])/(x2[i2]-x2[0])+interpolate1Cr(x2,y2,n2,&j,x3[i])*(x3[i]-x2[0])/(x2[i2]-x2[0]);
					y3[2*i+1]=y1[2*i+1]*(x2[i2]-x3[i])/(x2[i2]-x2[0])+interpolate1Ci(x2,y2,n2,&j,x3[i])*(x3[i]-x2[0])/(x2[i2]-x2[0]); }
				for(;i<n;i++) {y3[2*i]=y2[2*i2];y3[2*i+1]=y2[2*(i2++)+1];}}
			else {
				for(i=0;x3[i]<x2[0];i++) y3[i]=y1[i];
				for(j=-2;x3[i]<=x1[s1->n-1];i++)
					y3[i]=y1[i]*(x2[i2]-x3[i])/(x2[i2]-x2[0])+interpolate1(x2,y2,n2,&j,x3[i])*(x3[i]-x2[0])/(x2[i2]-x2[0]);
				for(;i<n;i++) y3[i]=y2[i2++]; }
			if(i2!=n2) er=999; }
	else if(!strcmp(fn,"log x")) {
		for(i=0;i<n;i++)	x3[i]=(x3[i]>0)?log10(x3[i]):0;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"ln x")) {
		for(i=0;i<n;i++)	x3[i]=(x3[i]>0)?log(x3[i]):0;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"exp x")) {
		for(i=0;i<n;i++)	x3[i]=exp(x3[i]);
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"10^x")) {
		for(i=0;i<n;i++)	x3[i]=pow(10,x3[i]);
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"sqrt x")) {
		for(i=0;i<n;i++)	x3[i]=(x3[i]>0)?sqrt(x3[i]):0;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"x+k") || !strcmp(fn,"k+x")) {
		for(i=0;i<n;i++)	x3[i]+=k;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"x-k")) {
		for(i=0;i<n;i++)	x3[i]-=k;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"k-x")) {
		for(i=0;i<n;i++)	x3[i]=k-x3[i];
		sort=1;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"x*k") || !strcmp(fn,"k*x")) {
		for(i=0;i<n;i++)	x3[i]*=k;
		if(k<0) sort=1;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"x/k")) {
		if(!k) {er=9;goto failure;}
		for(i=0;i<n;i++)	x3[i]/=k;
		if(k<0) sort=1;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"k/x")) {
		for(i=0;i<n;i++)	x3[i]=x3[i]?k/x3[i]:0;
		sort=1;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"x^k")) {
		for(i=0;i<n;i++)
			x3[i]=((x3[i]==0 && k<=0) || (x3[i]<0 && !isinteger(k)))?0:pow(x3[i],k);
		sort=1;
		copyV(y1,y3,cmplx?2*n:n); }
	else if(!strcmp(fn,"complex"))
		if(cmplx || (s2 && s2->cmplx)) {er=8;goto failure;}
		else {
			y3=allocV(2*n);
			if(!y3) {er=1;goto failure;}
			freeV(s3->y);
			s3->y=y3;
			s3->cmplx=1;
			for(i=0,j=-2;i<n;i++) {
				y3[2*i]=y1[i];
				y3[2*i+1]=s2?((x2[i]==x1[i])?y2[i]:interpolate1(x2,y2,n2,&j,x1[i])):0; }}
	else if(!strcmp(fn,"real"))
		if(cmplx) {
			y3=allocV(n);
			if(!y3) {er=1;goto failure;}
			freeV(s3->y);
			s3->y=y3;
			s3->cmplx=0;
			real(y1,y3,n); }
		else {er=12;goto failure;}
	else if(!strcmp(fn,"imag"))
		if(cmplx) {
			y3=allocV(n);
			if(!y3) {er=1;goto failure;}
			freeV(s3->y);
			s3->y=y3;
			s3->cmplx=0;
			imag(y1,y3,n); }
		else {er=12;goto failure;}
	else if(!strcmp(fn,"fourier"))
		if(cmplx) {
			if(!isevenspV(x1,n,0.001)) {er=15;goto failure;}
			f1=2*PI/((maxV(x1,n)-minV(x1,n))*n/(n-1));
			for(i=0;i<n;i++) x3[i]=k+i*f1;
			fourier(x1,y1,x3,y3,n,n,-1); }
		else {er=12;goto failure;}
	else if(!strcmp(fn,"invfourier"))
		if(cmplx) {
			if(!isevenspV(x1,n,0.001)) {er=15;goto failure;}
			f1=2*PI/((maxV(x1,n)-minV(x1,n))*n/(n-1));
			for(i=0;i<n;i++) x3[i]=k+i*f1;
			fourier(x1,y1,x3,y3,n,n,1); }
		else {er=12;goto failure;}
	else if(!strcmp(fn,"ftpower"))
		if(cmplx) {er=8;goto failure;}
		else {
			if(!isevenspV(x1,n,0.001)) {er=15;goto failure;}
			sz1=sz2=NULL;
			er=SpectMath2(s1,&sz2,num,0,"baseline","ends");
			if(er) goto failure;
			er=SpectMath(sz2,NULL,&sz1,"complex",0);
			if(er) {SpectFree(sz2);goto failure;}
			f1=(maxV(x1,n)-minV(x1,n))/n;
			if(n%2) er=SpectMath(sz1,NULL,&sz2,"fourier",PI/f1*(-1.0+1.0/n));
			else er=SpectMath(sz1,NULL,&sz2,"fourier",-PI/f1);
			if(er) {SpectFree(sz1);SpectFree(sz2);goto failure;}
			for(i=0;i<n;i++)	{
				x3[i]=sz2->x[i];
				y3[i]=CMPXmag(sz2->y[2*i],sz2->y[2*i+1]); }
			SpectFree(sz1);
			SpectFree(sz2); }
	else if(!strcmp(fn,"fft"))
		if(cmplx) {
			if(!isevenspV(x1,n,0.001)) {er=15;goto failure;}
			if(!is2ton(n)) {er=11;goto failure;}
			f1=2*PI/((maxV(x1,n)-minV(x1,n))*n/(n-1));	// Ęk
			for(i=0;i<n;i++) x3[i]=k+i*f1;
			fft(x1,y1,x3,y3,n,-1); }
		else {er=12;goto failure;}
	else if(!strcmp(fn,"invfft"))
		if(cmplx) {
			if(!isevenspV(x1,n,0.001)) {er=15;goto failure;}
			if(!is2ton(n)) {er=11;goto failure;}
			f1=2*PI/((maxV(x1,n)-minV(x1,n))*n/(n-1));
			for(i=0;i<n;i++) x3[i]=k+i*f1;
			fft(x1,y1,x3,y3,n,1); }
		else {er=12;goto failure;}
	else if(!strcmp(fn,"hankel"))
		if(cmplx) {er=8;goto failure;}
		else {
			if(!isevenspV(x1,n,0.001)) {er=15;goto failure;}
			f1=2*PI/((maxV(x1,n)-minV(x1,n))*n/(n-1));
			if(k<=0) k=10;
			for(i=0;i<n;i++) x3[i]=i*f1;
			hankel(x1,y1,x3,y3,n,n,k); }
	else if(!strcmp(fn,"noise"))
		if(cmplx) for(i=0;i<n;i++) {
				y3[2*i]=binomrand(5,0,k);
				y3[2*i+1]=binomrand(5,0,k); }
		else for(i=0;i<n;i++)	y3[i]=binomrand(5,0,k);
	else if(!strcmp(fn,"convolve")) {
		if(!s2) {er=6;goto failure;}
		if(cmplx || s2->cmplx) {er=8;goto failure;}
		if(!isevenspV(x1,n,0.001)) {er=15;goto failure;}
		f1=(x1[n-1]-x1[0])/(n-1);
		f2=(x2[n2-1]-x2[0])/(n2-1);
		j=floor(-x2[0]/f1+0.5);
		if(isevenspV(x2,n2,0.001) && fabs(f1/f2-1.0)<0.0001 && fabs(j+x2[0]/f1)<0.0001) {
			v1=x2;
			v2=y2;
			i2=n2; }
		else {
			i2=(x2[n2-1]-x2[0])/f1+1;
			v1=allocV(i2);
			if(!v1) {er=1;goto failure;}
			v2=allocV(i2);
			if(!v2) {freeV(v1);er=1;goto failure;}
			for(i=0;i<i2;i++) v1[i]=(i-j)*f1;
			convertxV(x2,y2,v1,v2,n2,i2);
			freeV(v1); }
		convolveV(y1,v2,y3,n,i2,n,j,y1[0],y1[n-1]);
		if(v2!=y2) freeV(v2); }
	else	{
		er=7;
		goto failure; }

	if(sort) {
		if(s3->cmplx) sortCV(s3->x,s3->y,s3->n);
		else sortV(s3->x,s3->y,s3->n); }
	*s3ptr=s3;
	return 0;

 failure:
  if(!*s3ptr) SpectFree(s3);
 	return er; }


int SpectMath2(sptr s,sptr *ansptr,float *num,int nn,char *fn,char *fn2) {
	sptr s2,s3,s4,ans;
	int i,j,n,i1,er;
	float x,f,num2[MAXARG];

	er=0;
	if(!ansptr) return 6;
	if(*ansptr) SpectFree(*ansptr);
	*ansptr=ans=s2=s3=s4=NULL;

	if(!fn || !fn2 || !num)
		return 6;
	else if(!strcmp(fn,"new")) {
		if(nn<3) return 26;
		if(nn>3) return 25;
		if(num[2]<=0) return 14;
		ans=SpectAlloc(NULL,NULL,fn2,"","");
		if(!ans) return 1;
		if(num[0]>num[1])	{f=num[0];num[0]=num[1];num[1]=f;}
		ans->n=(int) ((num[1]-num[0])/num[2])+1;
		ans->x=allocV(ans->n);
		ans->y=allocV(ans->n);
		if(!ans->x || !ans->y) {SpectFree(ans);return 1;}
		for(i=0;i<ans->n;i++)	ans->x[i]=i*num[2]+num[0];
		if(!strcmp(fn2,"1")) for(i=0;i<ans->n;i++)	ans->y[i]=1;
		else if(!strcmp(fn2,"x")) for(i=0;i<ans->n;i++)	ans->y[i]=ans->x[i];
		else {SpectFree(ans);return 13; }}
	else if(!s)
		return 6;
	else if(!strcmp(fn,"copy") && (nn<3 || num[2]<=0)) {
		if(nn<2) return 26;
		if(nn>3) return 25;
		if(num[0]>num[1])	{f=num[0];num[0]=num[1];num[1]=f;}
		i=locateV(s->x,num[0],s->n);
		if(i<0 || s->x[i]<num[0]) i++;
		j=locateV(s->x,num[1],s->n);
		if(i>j) return 11;
		n=s->n;
		s->n=j-i+1;
		s->x+=i;
		s->y+=s->cmplx?2*i:i;
		er=SpectMath(s,NULL,&ans,"copy",0);
		s->n=n;
		s->x-=i;
		s->y-=s->cmplx?2*i:i;
		if(er) return er; }
	else if(!strcmp(fn,"copy"))	{
		if(nn<3) return 26;
		if(nn>3) return 25;
		if(num[2]<=0) return 14;
		ans=SpectAlloc(NULL,NULL,s->desc,s->xunit,s->yunit);
		if(!ans) return 1;
		strcpy(ans->color,s->color);
		if(num[0]>num[1])	{f=num[0];num[0]=num[1];num[1]=f;}
		ans->n=(int) ((num[1]-num[0])/num[2])+1;
		ans->x=allocV(ans->n);
		ans->y=allocV(s->cmplx?2*ans->n:ans->n);
		if(!ans->x || !ans->y) {SpectFree(ans);return 1;}
		for(j=-2,i=0;i<ans->n;i++)	{
			ans->x[i]=i*num[2]+num[0];
			if(s->cmplx) {
				ans->y[2*i]=interpolate1Cr(s->x,s->y,s->n,&j,ans->x[i]);
				ans->y[2*i+1]=interpolate1Ci(s->x,s->y,s->n,&j,ans->x[i]); }
			else ans->y[i]=interpolate1(s->x,s->y,s->n,&j,ans->x[i]); }
		strncat(ans->desc," copy",STRCHAR-strlen(ans->desc)); }
	else if(!strcmp(fn,"zerofill"))	{
		if(nn<1) return 26;
		ans=SpectAlloc(NULL,NULL,s->desc,s->xunit,s->yunit);
		if(!ans) return 0;
		strcpy(ans->color,s->color);
		if(num[0]<=0)	ans->n=next2ton(s->n);
		else if(num[0]<1) {SpectFree(ans);return 14;}
		else ans->n=(int)(s->n*num[0]);
		ans->x=allocV(ans->n);
		ans->y=allocV(s->cmplx?2*ans->n:ans->n);
		if(!ans->x || !ans->y) {SpectFree(ans);return 1;}
		ans->cmplx=s->cmplx;
		f=(maxV(s->x,s->n)-minV(s->x,s->n))/(1.0*s->n-1);
		if(!strcmp(fn2,"right")) i1=0;
		else if(!strcmp(fn2,"left")) i1=ans->n-s->n;
		else if(!strcmp(fn2,"ends")) i1=(ans->n-s->n)/2;
		else {SpectFree(ans);return 13;}
		for(i=0;i<i1;i++)	{
			ans->x[i]=s->x[0]-(i1-i)*f;
			if(s->cmplx) ans->y[2*i]=ans->y[2*i+1]=0;
			else ans->y[i]=0;}
		for(i=i1;i<i1+s->n;i++)	{
			ans->x[i]=s->x[i-i1];
			if(s->cmplx) {ans->y[2*i]=s->y[2*(i-i1)];ans->y[2*i+1]=s->y[2*(i-i1)+1];}
			else ans->y[i]=s->y[i-i1];}
		for(i=i1+s->n;i<ans->n;i++)	{
			ans->x[i]=s->x[s->n-1]+(i-i1-s->n+1)*f;
			if(s->cmplx) ans->y[2*i]=ans->y[2*i+1]=0;
			else ans->y[i]=0;}
		strncat(ans->desc," zerofill",STRCHAR-strlen(ans->desc)); }
	else if(!strcmp(fn,"baseline"))	{
		if(nn>0) return 25;
		if(s->cmplx) return 8;
		er=SpectMath(s,NULL,&ans,"copy",0);
		if(er) return er;
		if(!strcmp(fn2,"ends"))	{
			n=s->n;
			num[0]=(s->y[n-1]-s->y[0])/(s->x[n-1]-s->x[0]);
			num[1]=(s->y[n-1]*s->x[0]-s->y[0]*s->x[n-1])/(s->x[0]-s->x[n-1]);
			for(i=0;i<n;i++) ans->y[i]=s->y[i]-(num[0]*s->x[i]+num[1]); }
		else if(!strcmp(fn2,"left")) {
			num[0]=s->y[0];
			for(i=0;i<s->n;i++)	ans->y[i]=s->y[i]-num[0]; }
		else if(!strcmp(fn2,"right"))	{
			num[0]=s->y[s->n-1];
			for(i=0;i<s->n;i++)	ans->y[i]=s->y[i]-num[0]; }
		else	{SpectFree(ans);return 13;}}
	else if(!strcmp(fn,"unbaseline"))	{
		if(nn<1) return 26;
		if(nn>2) return 25;
		if(s->cmplx) return 8;
		er=SpectMath(s,NULL,&ans,"copy",0);
		if(er) return er;
		if(!strcmp(fn2,"ends"))
			for(i=0;i<s->n;i++)	ans->y[i]=s->y[i]+(num[0]*s->x[i]+num[1]);
		else if(!strcmp(fn2,"left"))
			for(i=0;i<s->n;i++)	ans->y[i]=s->y[i]+num[0];
		else if(!strcmp(fn2,"right"))
			for(i=0;i<s->n;i++)	ans->y[i]=s->y[i]+num[0];
		else	{SpectFree(ans);return 13;}}
	else if(!strcmp(fn,"mask"))	{
		if(nn<1) return 26;
		if(s->cmplx) return 8;
		er=SpectMath(s,NULL,&ans,"copy",0);
		if(er) return er;
		if(!strcmp(fn2,"notch"))
			for(i=0;i<s->n;i++)	{
				x=fabs(s->x[i]);
				for(j=0;j<nn;j+=2)
					ans->y[i]*=1-exp(-(x-num[j])*(x-num[j])/(2*num[j+1]*num[j+1]));}
		else if(!strcmp(fn2,"gauss"))
			for(i=0;i<s->n;i++)
				ans->y[i]*=exp(-s->x[i]*s->x[i]/(2*num[0]*num[0]));
		else if(!strcmp(fn2,"highpass"))
			for(i=0;i<s->n;i++)	{
				x=fabs(s->x[i]);
				if(x<num[0]-num[1]/2)	ans->y[i]=0;
				else if(x<num[0]+num[1]/2) ans->y[i]*=(x-(num[0]-num[1]/2))/num[1];}
		else if(!strcmp(fn2,"lowpass"))
			for(i=0;i<s->n;i++)	{
				x=fabs(s->x[i]);
				if(x>num[0]+num[1]/2)	ans->y[i]=0;
				else if(x>num[0]-num[1]/2) ans->y[i]*=((num[0]+num[1]/2)-x)/num[1];}
		else	{SpectFree(ans);return 13;}}
	else if(!strcmp(fn,"filter"))	{
		if(s->cmplx) return 8;
		n=s->n;
		er=SpectMath2(s,&s2,num2,0,"baseline","ends");
		if(!er) er=SpectMath(s2,NULL,&s3,"complex",0);
		if(!er) er=SpectMath(s3,NULL,&s2,"fourier",-PI/(s->x[1]-s->x[0]));
		if(!er) er=SpectMath(s2,NULL,&s3,"real",0);
		if(!er) er=SpectMath(s2,NULL,&s4,"imag",0);
		if(!er) er=SpectMath2(s3,&s2,num,nn,"mask",fn2);
		if(!er) er=SpectMath2(s4,&s3,num,nn,"mask",fn2);
		if(!er) er=SpectMath(s2,s3,&s4,"complex",0);
		if(!er) er=SpectMath(s4,NULL,&s2,"invfourier",s->x[0]);
		if(!er) er=SpectMath(s2,NULL,&s3,"real",0);
		if(!er) er=SpectMath2(s3,&ans,num2,2,"unbaseline","ends");
		SpectFree(s2);
		SpectFree(s3);
		SpectFree(s4);
		if(er) return er; }
	else
		return 7;

	*ansptr=ans;
	return 0; }