#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "nrutil.h"

void dfpmin(double p[], int n, double gtol, int *iter, double *fret, 
	double (*func)(double []), void (*dfunc)(double [], double []));
int strcount(char *line, char c);
char *paulgetline(char *line, int nchars, FILE *fp);
double spectrumError( double orientation[] );
void errorGradient( double orientationIn[], double gradientOut[] );

double *spectrum;
double **localstd;
int *squash;
int nnnn, tflag=1;

#define MIN(A,B) ((A > B) ? (B) : (A))
#define PI 3.141592653589793
#define GTOL 1.0e-09

double theta[] = {PI/4, 0.05, 1.52, 0.05, PI/4, 1.52, 0.05, 0.5, 1.0, 1.52};
double phi[] = {0.05, 0.5, 0.5, 1.0, PI/4, 1.0, 1.52, 1.52, 1.52, 1.52};

int main(int argc, char *argv[]) {
  FILE *fp1, *fp2;
  int q, i, j, ni, n1, n2, nstd=0, *arglist;
  int nflag=0, iflag=0, oflag=0, no, o=0, hl=2, iter, nxcl=0, yflag=0;
  double *x, *xstd, x1=1200.0, x2=2200.0, hold, lowxcl[20], highxcl[20];
  double **spectra, **stds, **stdinterp, **p, pp[4], *thickness, *fret, frettemp, intlow=3000.0, inthigh=3750.0;
  char file1[256], sfile[256], ofile[256], yfile[256], *line, *linestart;
  double *w, **v, **a, *y, grate[3], **synthetics;
  int *yes;
  
  if (argc == 1) {
    printf("Usage: pauldeconvolve [-i] [-t] [-n] [-s stdfilename] [-o outputfilename]\n");
	printf("        [-hl headerlines] [-1 lowwavenumber] [-2 highwavenumber]\n");
	printf("		[-x low1 high1 [-x low2 high2 ] ... ] [-y synthfile] [-int low high] [infile] \n");
	printf("-i puts you in interactive mode, giving you a chance to set all other parameters as prompted.\n");
	printf("-t indicates fixed thickness, otherwise thickness is varied during the fitting.\n");
	printf("-n indicates spectra are normalized and must be denormalized before fitting.\n");
	printf("-s specifies the standards file [no header lines!], (wavenumber, a, b, c) in csv format at 1 mm thick.\n");
	printf("-o specifies the name of the output file. By default .decon is appended to the input name.\n");
	printf("-hl specfies how many header lines in input file; these are copied to output file.\nThe line after the headerlines should give the thicknesses in mm.\n");
	printf("The spectrum is fitted between lowwavenumber and highwavenumber (defaults 1200-2200).\n");
	printf("-x indicates a region between low and high limits to be ignored in the fitting.\n");
	printf("-y gives the name of the synthesized (a,b,c) output file and turns on synth-ing.\n");
	printf("-int gives the wavenumber range to be integrated in the synthesized spectra (default 3000-3750).\n");
	exit(0);
  }
  arglist = (int *) malloc((unsigned) 300*sizeof(int));
  line = (char *) malloc((unsigned) 5000*sizeof(char));
  for (i=1;i<argc;i++) arglist[i] = 1;
  for (i=1;i<argc;i++) {
    if (!strcmp(argv[i], "-i")) {
	  iflag = 1; arglist[i] = 0;
	} else if (!strcmp(argv[i], "-n")) {
	  nflag = 1; arglist[i] = 0;
	} else if (!strcmp(argv[i], "-o")) {
	  oflag = 1; strcpy(ofile, argv[i+1]);
	  arglist[i++] = 0; arglist[i] = 0;
	} else if (!strcmp(argv[i], "-s")) {
	  strcpy(sfile, argv[i+1]);
	  arglist[i++] = 0; arglist[i] = 0;
	} else if (!strcmp(argv[i], "-hl")) {
      hl = atof(argv[i+1]);
	  arglist[i++] = 0; arglist[i] = 0;
	} else if (!strcmp(argv[i],"-1")) {
	  x1 = atof(argv[i+1]);
	  arglist[i++] = 0;
	  arglist[i] = 0;
    } else if (!strcmp(argv[i],"-2")) {
	  x2 = atof(argv[i+1]);
	  arglist[i++] = 0;
	  arglist[i] = 0;
	} else if (!strcmp(argv[i],"-x")) {
	  lowxcl[nxcl] = atof(argv[i+1]);
	  highxcl[nxcl++] = atof(argv[i+2]);
	  arglist[i++] = 0; arglist[i++] = 0; arglist[i] = 0;
	} else if(!strcmp(argv[i], "-t")) {
	  tflag = 0;
	  arglist[i] = 0;
	} else if (!strcmp(argv[i], "-y")) {
	  strcpy(yfile, argv[i+1]); yflag = 1;
	  arglist[i++] = 0; arglist[i] = 0;
	} else if (!strcmp(argv[i], "-int")) {
	  intlow = atof(argv[i+1]);
	  inthigh = atof(argv[i+2]);
	  arglist[i++] = 0; arglist[i++] = 0; arglist[i] = 0;
	}
  }
  for (q=1,o=0;q<argc;q++) {
    if (arglist[q]) {
	  strcpy(file1, argv[q]);
	  o++;
	}
  }
  if (o==0) {
    printf("No input file found on command line, going interactive...\n");
	iflag = 1;
  }
  
  if (iflag) {
    printf("Welcome to pauldeconvolve interactive mode!\n");
	printf("Step 1a. Enter standards filename [no header lines!] (or 'x' to use standards.csv): ");
	scanf("%s", sfile);
	if (!strcmp(sfile, "x")) strcpy(sfile, "standards.csv");
	printf("Step 1b. Enter output filename (or 'x' for <inputfile>.decon): ");
	scanf("%s", ofile);
	if (strcmp(ofile, "x")) oflag = 1;
	else oflag = 0;
	printf("Step 2a. Disallow thickness ratio variations? (0 for fixed thickness): ");
	scanf("%d", &tflag);
	printf("Step 2b. Are unknown spectra already normalized? (0 for as-measured, 1 for normalized): ");
	scanf("%d", &nflag);
	printf("Step 3a. Enter low wavenumber limit of fitted range (enter 0 for default: %f): ", x1);
	scanf("%lf", &hold); if (hold != 0.0) x1 = hold;
	printf("Step 3b. Enter low wavenumber limit of fitted range (enter 0 for default: %f): ", x2);
	scanf("%lf", &hold); if (hold != 0.0) x2 = hold;
    printf("Step 4. How many regions to exclude (now there are %d; 0 for no more): ", nxcl);
	scanf("%d", &q);
	for (i=nxcl;i<nxcl+q;i++) {
	  printf("Step 4-%da. Low wavenumber limit of excluded region %d: ", i+1, i+1);
	  scanf("%lf", &lowxcl[i]);
	  printf("Step 4-%db. High wavenumber limit of excluded region %d: ", i+1, i+1);
	  scanf("%lf", &highxcl[i]);
    }
	nxcl += q;
	printf("Step 5. Input filename: ");
	scanf("%s", file1);
	printf("Step 6. Synthesize principal spectra? (1 for yes, 0 for no): ");
	scanf("%d", &yflag);
	if (yflag) {
	  printf("Step 6a. Synthetics output filename: ");
	  scanf("%s", yfile);
	  printf("Step 6b. Low wavenumber limit of region to integrate in synthetics (0 for 3000): ");
	  scanf("%lf", &intlow); if (intlow == 0.0) intlow = 3000.0;
	  printf("Step 6b. High wavenumber limit of region to integrate in synthetics (0 for 3750): ");
	  scanf("%lf", &intlow); if (intlow == 0.0) intlow = 3750.0;
	}
  } 

  if ((fp1 = fopen(sfile, "r")) == NULL) {
	printf("Could not open standards file %s\n", sfile);
	exit(-1);
  }
  printf("Reading standards file %s\n", sfile);
  while ( paulgetline(line, 5000, fp1) != NULL) nstd++;
  rewind(fp1);
  xstd = (double *) malloc((unsigned) nstd*sizeof(double));
  stds = (double **) dmatrix(0,2,0,nstd);
  for (i=0;i<nstd;i++) {
    paulgetline(line, 5000, fp1);
	sscanf(line, "%lf,%lf,%lf,%lf", &xstd[i], &stds[0][i], &stds[1][i], &stds[2][i]);
  }
  fclose(fp1);

  if ((fp1 = fopen(file1, "r")) == NULL) {
	printf("Could not open file %s\n", file1);
	exit(-1);
  }
  printf("Reading file %s\n", file1);
  if (!oflag) {
	strcpy(ofile, file1);
	ofile[strlen(ofile)-3] = '\0';
	strcat(ofile, "decon.csv");
  }
  if ((fp2 = fopen(ofile, "w")) == NULL) {
    printf("Could not write file %s\n", ofile);
	exit(-1);
  }
  printf("Writing file %s\n", ofile);
  ni = -(hl+1); 
  while ( paulgetline(line, 5000, fp1) != NULL) {
	if (ni < -1) fprintf(fp2, "%s\n", line);
	if (ni == 0) no = strcount(line, ',');
	ni++;
  }  
  rewind(fp1);
  x = (double *) malloc((unsigned) ni*sizeof(double));
  p = (double **) dmatrix(0,no,1,3);
  fret = (double *) malloc((unsigned) no*sizeof(double));
  spectra = (double **) dmatrix(0,no,0,ni);
  thickness = (double *) malloc((unsigned) no*sizeof(double));
  w = (double *) malloc((unsigned) 4*sizeof(double));
  synthetics = (double **) dmatrix(0,ni,1,3);
  v = (double **) dmatrix(1,3,1,3);
  a = (double **) dmatrix(1,no,1,3);
  yes = (int *) malloc((unsigned) (no+1)*sizeof(int));
  y = (double *) malloc((unsigned) (no+1)*sizeof(double));
  for (i=1; i<=no; i++) yes[i] = 1;
  for (i=-(hl+1);i<ni;i++) {
	paulgetline(line, 5000, fp1);
	if (i==-1) {
	  linestart = line;
	  for (j=0;j<no;j++) {
		line = strchr(line, ',') + 1;
		sscanf(line, "%lf", &thickness[j]);
	  }
	  line = linestart;
	} else if (i >= 0) {
	  linestart = line;
	  sscanf(line, "%lf", &x[i]);
	  for (j=0;j<no;j++) {
		line = strchr(line, ',') + 1;
		sscanf(line, "%lf", &spectra[j][i]);
	  }
	  line = linestart;
	}
  }
  fclose(fp1);
  n1 = 0;
  while (x[n1] < x1) n1++;
  n2 = 0;
  while (x[n2] < x2) n2++;
  nnnn = n2 - n1;
  stdinterp = (double **) dmatrix(0,2,0,ni);
  localstd = (double **) dmatrix(0,2,0,nnnn);
  spectrum = (double *) dvector(0,nnnn);
  squash = (int *) ivector(0,nnnn);
  for (i=0;i<nnnn;i++) squash[i] = 0;
  for (q=0;q<nxcl;q++) {
    for (i=0;i<nnnn;i++) if (lowxcl[q] <= x[n1+i] && highxcl[q] >= x[n1+i]) squash[i] = 1;
  }
	  
  /* First, interpolate standards to sample x points */
  for (i=n1;i<n2;i++) {
    int sn = 0;
    while (xstd[sn] < x[i]) sn++;
    for (j=0; j<3; j++) {
	  stdinterp[j][i] = (stds[j][sn+1]*(x[i]-xstd[sn])+stds[j][sn]*(xstd[sn+1]-x[i]))/(xstd[sn+1]-xstd[sn]);
	}
  }
	
  /* convert to transmission; if spectra are normalized, denormalize them */
  for (j=0;j<no;j++) {
    for (i=0;i<ni;i++) {
	  if (nflag) spectra[j][i] *= thickness[j];
      spectra[j][i] = pow(10.0, -spectra[j][i]);
    }
  }
  
  for (j=0;j<no;j++) {  /* main loop over spectra */
	/* denormalize standards */
	for (q=0;q<3;q++) 
	  for (i=0;i<nnnn;i++) localstd[q][i] = stdinterp[q][i+n1]*thickness[j];
	
	/* convert standards to transmission, read in local spectrum */
	for (i=0;i<nnnn;i++) {
	  spectrum[i] = spectra[j][i+n1];
	  for (q=0;q<3;q++) localstd[q][i] = pow(10.0, -localstd[q][i]);
	}
	
	/* differentiate spectrum */
	for (i=0;i<nnnn-9;i++) spectrum[i] = (spectrum[i+9]-spectrum[i])/9.0;
	
	/* initial guess */
	for (fret[j] = 10000.0,i=0;i<(tflag? 30 : 10);i++) {
	  pp[1] = theta[i%10];
	  pp[2] = phi[i%10];
	  if (tflag) pp[3] = (i<10 ? 0.9: (i<20 ? 1.0 : 1.1));
	  else pp[3] = 1.0;
	  dfpmin(pp, 3, GTOL, &iter, &frettemp, spectrumError, errorGradient);
	  if (frettemp < fret[j]) {
	    p[j][1] = asin(sin(acos(cos(pp[1]))));
	    p[j][2] = asin(sin(acos(cos(pp[2]))));
	    p[j][3] = pp[3];
	    fret[j] = frettemp;
	  }
	}
	printf("Spectrum %d: Converged to %g in %d iterations: Th = %g, Ph = %g, Tr = %f\n", j, fret[j], iter, (180.0/PI)*p[j][1], (180.0/PI)*p[j][2], p[j][3]);
  }

  fprintf(fp2, "Theta");
  for (j=0;j<no;j++) fprintf(fp2, ",%f", (180.0/PI)*p[j][1]);
  fprintf(fp2, "\nPhi");
  for (j=0;j<no;j++) fprintf(fp2, ",%f", (180.0/PI)*p[j][2]);
  fprintf(fp2, "\nThicknessRatio");
  for (j=0;j<no;j++) fprintf(fp2, ",%f", p[j][3]);
  fprintf(fp2, "\nError");
  for (j=0;j<no;j++) fprintf(fp2, ",%f", fret[j]);
  fclose(fp2);
  
  if (yflag) {			
    if ((fp2 = fopen(yfile, "w")) == NULL) {
	  printf("Could not open synthesized principal spectra output file.\n");
	  exit(-1);
	}

    /* make orientation matrix */
	for (i=1;i<=no;i++) {
	  a[i][1] = cos(p[i-1][1])*cos(p[i-1][1])*sin(p[i-1][2])*sin(p[i-1][2]);
	  a[i][2] = sin(p[i-1][1])*sin(p[i-1][1])*sin(p[i-1][2])*sin(p[i-1][2]);
	  a[i][3] = cos(p[i-1][2])*cos(p[i-1][2]);
	}
	/* svd */
	dsvdcmp(a, no, 3, w, v);
	for (j=1;j<=3;j++) if (w[j] < 1.0e-08) w[j] = 0;
	
	/* renormalize to 0.1 mm and backsub at each wavenumber */
	for (i=0;i<ni;i++) {
	  for (j=0;j<no;j++) y[j+1] = pow(spectra[j][i], 0.01/(thickness[j]));
	  dsvbksb(a, w, v, no, 3, y, synthetics[i]);
	}			

    /* convert to absorbance and renormalize to 1 cm*/
	for (i=0;i<ni;i++)
	  for (j=1;j<=3;j++) synthetics[i][j] = -1000.0*log10(synthetics[i][j]);
	  
    /* Integrate the synthetics */
	n1 = 0;
    while (x[n1] < intlow) n1++;
    n2 = 0;
    while (x[n2] < inthigh) n2++;
	for (j=1;j<=3;j++) {
	  for (grate[j]=0.0,i=n1;i<n2;i++) grate[j] += synthetics[i][j];
	  grate[j] *= (x[n2]-x[n1])/(n2-n1);
	}
	
	/* output */
	fprintf(fp2, "orientation:,a,b,c,total\n");
	fprintf(fp2, "integral:,%g,%g,%g,%g\n", grate[1], grate[2], grate[3], grate[1]+grate[2]+grate[3]);
	fprintf(fp2, "water:,%g,%g,%g,%g\n", 0.188*grate[1], 0.188*grate[2], 0.188*grate[3], 0.188*(grate[1]+grate[2]+grate[3]));
	for (i=0;i<ni;i++) {
	  fprintf(fp2, "%f,%g,%g,%g\n", x[i], synthetics[i][1], synthetics[i][2], synthetics[i][3]);
	}
	fclose(fp2);
  }
  return 0;
}


#define ITMAX 500
#define EPS 3.0e-8
#define TOLX (4*EPS)
#define STPMX 100.0

#define FREEALL free_dvector(xi,1,n);free_dvector(pnew,1,n); \
free_dmatrix(hessin,1,n,1,n);free_dvector(hdg,1,n);free_dvector(g,1,n); \
free_dvector(dg,1,n);

void dfpmin(double p[], int n, double gtol, int *iter, double *fret,
	double(*func)(double []), void (*dfunc)(double [], double []))
{
	void lnsrch(int n, double xold[], double fold, double g[], double p[], double x[],
		 double *f, double stpmax, int *check, double (*func)(double []));
	int check,i,its,j;
	double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test;
	double *dg,*g,*hdg,**hessin,*pnew,*xi;

	dg=dvector(1,n);
	g=dvector(1,n);
	hdg=dvector(1,n);
	hessin=dmatrix(1,n,1,n);
	pnew=dvector(1,n);
	xi=dvector(1,n);
	fp=(*func)(p);
	(*dfunc)(p,g);
	for (i=1;i<=n;i++) {
		for (j=1;j<=n;j++) hessin[i][j]=0.0;
		hessin[i][i]=1.0;
		xi[i] = -g[i];
		sum += p[i]*p[i];
	}
	stpmax=STPMX*FMAX(sqrt(sum),(double)n);
	for (its=1;its<=ITMAX;its++) {
		*iter=its;
		lnsrch(n,p,fp,g,xi,pnew,fret,stpmax,&check,func);
		fp = *fret;
		for (i=1;i<=n;i++) {
			xi[i]=pnew[i]-p[i];
			p[i]=pnew[i];
		}
		test=0.0;
		for (i=1;i<=n;i++) {
			temp=fabs(xi[i])/FMAX(fabs(p[i]),1.0);
			if (temp > test) test=temp;
		}
		if (test < TOLX) {
			FREEALL
			return;
		}
		for (i=1;i<=n;i++) dg[i]=g[i];
		(*dfunc)(p,g);
		test=0.0;
		den=FMAX(*fret,1.0);
		for (i=1;i<=n;i++) {
			temp=fabs(g[i])*FMAX(fabs(p[i]),1.0)/den;
			if (temp > test) test=temp;
		}
		if (test < gtol) {
			FREEALL
			return;
		}
		for (i=1;i<=n;i++) dg[i]=g[i]-dg[i];
		for (i=1;i<=n;i++) {
			hdg[i]=0.0;
			for (j=1;j<=n;j++) hdg[i] += hessin[i][j]*dg[j];
		}
		fac=fae=sumdg=sumxi=0.0;
		for (i=1;i<=n;i++) {
			fac += dg[i]*xi[i];
			fae += dg[i]*hdg[i];
			sumdg += SQR(dg[i]);
			sumxi += SQR(xi[i]);
		}
		if (fac*fac > EPS*sumdg*sumxi) {
			fac=1.0/fac;
			fad=1.0/fae;
			for (i=1;i<=n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
			for (i=1;i<=n;i++) {
				for (j=1;j<=n;j++) {
					hessin[i][j] += fac*xi[i]*xi[j]
					-fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
				}
			}
		}
		for (i=1;i<=n;i++) {
			xi[i]=0.0;
			for (j=1;j<=n;j++) xi[i] -= hessin[i][j]*g[j];
		}
	}
	nrerror("too many iterations in dfpmin");
	FREEALL
}
#undef ITMAX
#undef EPS
#undef TOLX
#undef STPMX
#undef FREEALL
#define ALF 1.0e-4
#define TOLX 1.0e-7

void lnsrch(int n, double xold[], double fold, double g[], double p[], double x[],
	double *f, double stpmax, int *check, double (*func)(double []))
{
	int i;
	double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,
		test,tmplam;

	*check=0;
	for (sum=0.0,i=1;i<=n;i++) sum += p[i]*p[i];
	sum=sqrt(sum);
	if (sum > stpmax)
		for (i=1;i<=n;i++) p[i] *= stpmax/sum;
	for (slope=0.0,i=1;i<=n;i++)
		slope += g[i]*p[i];
	test=0.0;
	for (i=1;i<=n;i++) {
		temp=fabs(p[i])/FMAX(fabs(xold[i]),1.0);
		if (temp > test) test=temp;
	}
	alamin=TOLX/test;
	alam=1.0;
	for (;;) {
		for (i=1;i<=n;i++) x[i]=xold[i]+alam*p[i];
		*f=(*func)(x);
		if (alam < alamin) {
			for (i=1;i<=n;i++) x[i]=xold[i];
			*check=1;
			return;
		} else if (*f <= fold+ALF*alam*slope) return;
		else {
			if (alam == 1.0)
				tmplam = -slope/(2.0*(*f-fold-slope));
			else {
				rhs1 = *f-fold-alam*slope;
				rhs2=f2-fold2-alam2*slope;
				a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
				b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
				if (a == 0.0) tmplam = -slope/(2.0*b);
				else {
					disc=b*b-3.0*a*slope;
					/*if (disc<0.0) nrerror("Roundoff problem in lnsrch.");
					  else */tmplam=(-b+sqrt(disc))/(3.0*a);
				}
				if (tmplam>0.5*alam)
					tmplam=0.5*alam;
			}
		}
		alam2=alam;
		f2 = *f;
		fold2=fold;
		alam=FMAX(tmplam,0.1*alam);
	}
}
#undef ALF
#undef TOLX

int strcount(char *line, char c) {
  int count=0;
  char *marker = line;
  while ((marker = strchr(marker, c)) != '\0') {count++; marker++;}
  return count;
}

char *paulgetline(char *line, int nchars, FILE *fp) {
  char c, *startline=line;
  int n=0;

  while ((c = fgetc(fp)) != EOF && c != 10 && c !=13 && n++ <= nchars) *line++ = c;
  if (c == 13) { 
	if ((c = fgetc(fp)) != EOF && c != 10) {
	  ungetc(c, fp);
	}
  }
  *line = '\0';
  if (line == startline) return 0x0;
  else {
    line = startline;
    return line;
  }
}

double spectrumError( double orientation[] ) {
  double theta = orientation[1];
  double phi = orientation[2];
  double tratio = orientation[3];
  double *test, errorVal, errorInc, a, b, c;
  int i;
  
  test = (double *) malloc((unsigned) nnnn*sizeof(double));
  /* generate test spectrum */
  a = cos(theta)*sin(phi); a *= a;
  b = sin(theta)*sin(phi); b *= b;
  c = cos(phi); c *= c;
  for (i=0;i<nnnn;i++)
    test[i] = pow(localstd[0][i],tratio)*a+pow(localstd[1][i],tratio)*b+pow(localstd[2][i],tratio)*c;
  /* differentiate it */
  for (i=0;i<nnnn-9;i++)
  	test[i] = (test[i+9]-test[i])/9.0;
  /* get the error */
  for (errorVal=0,i=0;i<nnnn-9;i++) {
    if (!squash[i]) {
      errorInc = spectrum[i] - test[i];
  	  errorInc *= errorInc;
	  errorVal += errorInc;
	}
  }
	
  free(test);
  return errorVal;
}

void errorGradient( double orientationIn[], double gradientOut[] ) {
  double theta = orientationIn[1];
  double phi = orientationIn[2];
  double tratio = orientationIn[3];
  double *test, errorInc, a, b, c, ab, ac, bc, cc;
  int i;
  
  test = (double *) malloc((unsigned) nnnn*sizeof(double));
  /* generate test spectrum */
  a = cos(theta)*sin(phi); a *= a;
  b = sin(theta)*sin(phi); b *= b;
  c = cos(phi); c *= c;
  ab = cos(theta)*sin(theta)*sin(phi)*sin(phi);
  ac = cos(theta)*cos(theta)*sin(phi)*cos(phi);
  bc = sin(theta)*sin(theta)*sin(phi)*cos(phi);
  cc = sin(phi)*cos(phi);
  for (i=0;i<nnnn;i++)
    test[i] = pow(localstd[0][i],tratio)*a+pow(localstd[1][i],tratio)*b+pow(localstd[2][i],tratio)*c;
  /* differentiate it */
  for (i=0;i<nnnn-9;i++)
  	test[i] = (test[i+9]-test[i])/9.0;
  /* get the error */
  for (gradientOut[1]=0.0,gradientOut[2]=0.0,gradientOut[3]=0.0,i=0;i<nnnn-9;i++) {
    errorInc = spectrum[i] - test[i];
    gradientOut[1] += errorInc*2.0*((pow(localstd[1][i],tratio)-pow(localstd[0][i],tratio))-(pow(localstd[1][i+9],tratio)-pow(localstd[0][i+9],tratio)))*ab;
	gradientOut[2] += errorInc*2.0*(ac*(pow(localstd[0][i],tratio)-pow(localstd[0][i+9],tratio))+\
									bc*(pow(localstd[1][i],tratio)-pow(localstd[1][i+9],tratio))-\
									cc*(pow(localstd[2][i],tratio)-pow(localstd[2][i+9],tratio)));
	if (tflag) gradientOut[3] += errorInc*(a*(pow(localstd[0][i],tratio)*log(localstd[0][i])-pow(localstd[0][i+9],tratio)*log(localstd[0][i+9]))+\
								b*(pow(localstd[1][i],tratio)*log(localstd[1][i])-pow(localstd[1][i+9],tratio)*log(localstd[1][i+9]))+\
								c*(pow(localstd[2][i],tratio)*log(localstd[2][i])-pow(localstd[2][i+9],tratio)*log(localstd[2][i+9])));
	else gradientOut[3] = 0.0;
  }
  gradientOut[1] *= 2.0/9.0; gradientOut[2] *= 2.0/9.0; gradientOut[3] *= 2.0/9.0;	
  free(test);
  return;
}