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

void spline(double x[], double y[], int n, double yp1, double ypn, double y2[]);
void splint(double xa[], double ya[], double y2a[], int n, double x, double *y);
char *paulgetline(char *line, int nchars, FILE *fp);

#define MIN(A,B) ((A > B) ? (B) : (A))

int main(int argc, char *argv[]) {
  FILE *fp1, *fp2;
  int q, i, j, k, n, n1, n2, a[101], z, *arglist, nanchors=25, z1, z2, pz1, pz2, nf=0, nz=0, xed, Zflag=0, epoxypeakn;
  int npz=0, bgflag=0, ng=0, pflag=0, pdegree=5, gflag, iflag=0, zarg, oflag=0, oiflag=0, no, o=0, nx=0, ne;
  double *x, *spec, *back, x1=1100.0, x2=3970.0, delta, hold, xray[101], y[101], yy[101], max, swap, x00;
  double zap1[20], zap2[20], pzap1[20], pzap2[20], xr1[20], xr2[20], f[20], g[20], **oback, **iback;
  double *epoxyx, *epoxyy, epoxypeakx=2935.126, epoxyscale;
  char file1[256], file2[256], ofile[256], ifile[256], line[100], c[20], allfiles[5000], epoxyfile[256];
  
  if (argc == 1) {
    printf("Usage: pbg [-i] [-bgs] [-o omnibusfilename] [-oi omnibusinputfile] [-p degree]\n");
	printf("         [-1 lowwavenumber] [-2 highwavenumber] [-a #anchors] [-z low1 high1 [-z low2 high2] ... ]\n");
	printf("		 [-pz low1 high1 [-pz low2 high2] ... ] [-f fixedanchor1 [-f fixedanchor2] ... ]\n");
	printf("		 [-x low 1 high1 [-x low2 high2] ... ] [-g linearpoint1 [-g linearpoint2] ... ]\n");
	printf("		 [-Z epoxyfilename] [-Zp peaklocation] infile1 [infile2] [infile3]...\n");
	printf("-i puts you in interactive mode, giving you a chance to set all other parameters as prompted.\n");
	printf("-bgs means output the subtracted spectrum, otherwise the background itself is output.\n");
	printf("-o means assemble all output spectra into one omnibus file, in addition to the separate files.\n");
	printf("-oi means echo all input spectra into a single file, for your convenience.\n");
	printf("-p means polynomial fit, otherwise cubic spline is used.\n");
	printf("The background is fitted between lowwavenumber and highwavenumber, using #anchors anchor points.\n");
	printf("The -z parameter indicates a subregion of organic peaks to be copied into background after fitting.\n");
	printf("The -pz parameter indicates really troublesome regions to be replaced with lines before fitting.\n");
	printf("The -f points are wavenumbers of manually fixed anchors.\n");
	printf("The -x ranges are forbidden areas - the program will not pick anchors in these regions.\n");
	printf("Wavenumber ranges containing a g-point will be fitted with lines, not splines.\n");
	printf("-Z indicates that epoxy contamination spectrum epoxyfilename is to be subtracted.\n");
	printf("If -Zp is given then the epoxy intensity is determined at peaklocation, otherwise at 3925.126.\n");
	printf("Defaults: -bg, lowwavenumber = 1100, highwavenumber = 3970, #anchors = 25, no zap, no fixed anchors\n");
	printf("Note: Please ensure that all spectra run together in a single execution are sampled at exactly the same wavenumbers!\n");
	exit;
  }
  arglist = (int *) malloc((unsigned) 300*sizeof(int));
  for (i=1;i<argc;i++) arglist[i] = 1;
  for (i=1;i<argc;i++) {
    if (!strcmp(argv[i], "-bgs")) {
	  bgflag = 1; arglist[i] = 0;
	} else if (!strcmp(argv[i], "-i")) {
	  iflag = 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], "-oi")) {
	  oiflag = 1; strcpy(ifile, argv[i+1]);
	  arglist[i++] = 0; arglist[i] = 0;
	} else if (!strcmp(argv[i],"-p")) {
	  pflag = 1;
	  pdegree = atoi(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], "-a")) {
	  nanchors = atoi(argv[i+1]);
	  arglist[i++] = 0;
	  arglist[i] = 0;
	} else if(!strcmp(argv[i], "-z")) {
	  zap1[nz] = atof(argv[i+1]);
	  zap2[nz++] = atof(argv[i+2]);
	  arglist[i++] = 0; arglist[i++] = 0; arglist[i] = 0;
	} else if(!strcmp(argv[i], "-pz")) {
	  pzap1[npz] = atof(argv[i+1]);
	  pzap2[npz++] = atof(argv[i+2]);
	  arglist[i++] = 0; arglist[i++] = 0; arglist[i] = 0;
	} else if(!strcmp(argv[i], "-x")) {
	  xr1[nx] = atof(argv[i+1]);
	  xr2[nx++] = atof(argv[i+2]);
	  arglist[i++] = 0; arglist[i++] = 0; arglist[i] = 0;
	} else if(!strcmp(argv[i], "-f")) {
	  f[nf++] = atof(argv[i+1]);
	  arglist[i++] = 0; arglist[i] = 0;
	} else if(!strcmp(argv[i], "-g")) {
	  g[ng++] = atof(argv[i+1]);
	  arglist[i++] = 0; arglist[i] = 0;
	} else if(!strcmp(argv[i], "-Z")) {
	  Zflag=1; strcpy(epoxyfile, argv[i+1]);
	  arglist[i++] = 0; arglist[i] = 0;
	} else if(!strcmp(argv[i], "-Zp")) {
	  epoxypeakx = atof(argv[i+1]);
	  arglist[i++] = 0; arglist[i] = 0;
	}
  }
  
  if (iflag) {
    printf("Welcome to paulbackground2 interactive mode!\n");
	printf("Step 1a. Enter 'y' to output the subtracted spectrum (default: output the background): ");
	scanf("%s", &c);
	if (*c=='y') bgflag = 1;
	printf("Step 1b. Enter 'y' to output an omnibus file (default: only the separate output files): ");
	scanf("%s", &c);
	if (*c == 'y') {
	  oflag = 1;
	  printf("Step 1bb. Omnibus output filename: ");
	  scanf("%s", ofile);
	}
	printf("Step 1c. Enter 'y' to assemble inputs into an omnibus file (default: don't do this): ");
	scanf("%s", &c);
	if (*c == 'y') {
	  oiflag = 1;
	  printf("Step 1cc. Omnibus input spectra filename: ");
	  scanf("%s", ifile);
	}
	printf("Step 2. Enter 'y' for polynomial fit (default: cubic spline fit): ");
	scanf("%s", &c);
	if (*c=='y') {
	  pflag = 1;
	  printf("Step 2a. Degree of polynomial fit (enter 0 for default: %d): ", pdegree);
	  scanf("%d", &n); if (n != 0) pdegree = n;
	}
	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. Enter maximum number of anchors (enter 0 for default %d): ", nanchors);
	scanf("%d", &n); if (n != 0) nanchors = n;
	printf("Step 5. How many fixed anchor points do you want to add (now there are %d; enter 0 for no more): ", nf);
	scanf("%d", &n);
	for (i=0;i<n;i++) {
	  printf("Step 5-%d. Wavenumber of new fixed anchor number %d: ", i+1, i+1);
	  scanf("%lf", &f[nf+i]);
	  if (f[nf+i] < x1 || f[nf+i] > x2) {
	    printf("Outside fitted range, rejected. Try again.\n");
		i--;
	  }
	}
	nf += n;
	if (!pflag) {
      printf("Step 6. How many linear (non-spline) segments do you want to define (now there are %d; enter 0 for no more): ", ng);
	  scanf("%d", &n);
	  for (i=0;i<n;i++) {
	    printf("Step 6-%d. Wavenumber within new linear segment %d: ", i+1, i+1);
	    scanf("%lf", &g[nf+i]);
	    if (g[nf+i] < x1 || g[nf+i] > x2) {
	      printf("Outside fitted range, rejected. Try again.\n");
		  i--;
	    }
	  }
	  ng += n;
	}
	printf("Step 7. How many regions to copy peaks into background (now there are %d; enter 0 for no more): ", nz);
	scanf("%d", &n);
	for (i=0;i<n;i++) {
	  printf("Step 7-%da. Low wavenumber limit of new zap region %d: ", i+1, i+1);
	  scanf("%lf", &zap1[nz+i]);
	  printf("Step 7-%db. High wavenumber limit of new zap region %d: ", i+1, i+1);
	  scanf("%lf", &zap2[nz+i]);
	  if (zap1[nz+i] < x1 || zap1[nz+i] > x2 || zap2[nz+i] < x1 || zap2[nz+i] > x2) {
	    printf("One of those limits is outside fitted range, rejected. Try again.\n");
		i--;
	  }
	}
	nz += n;
	printf("Step 8. How many regions to replace with lines before fitting (now there are %d; enter 0 for no more): ", npz);
	scanf("%d", &n);
	for (i=0;i<n;i++) {
	  printf("Step 8-%da. Low wavenumber limit of new pre-zap region %d: ", i+1, i+1);
	  scanf("%lf", &pzap1[npz+i]);
	  printf("Step 8-%db. High wavenumber limit of new pre-zap region %d: ", i+1, i+1);
	  scanf("%lf", &pzap2[npz+i]);
	  if (pzap1[npz+i] < x1 || pzap1[npz+i] > x2 || pzap2[npz+i] < x1 || pzap2[npz+i] > x2) {
	    printf("One of those limits is outside fitted range, rejected. Try again.\n");
		i--;
	  }
	}
	npz += n;
	printf("Step 9. How many regions to forbid from anchoring (now there are %d; enter 0 for no more): ", nx);
	scanf("%d", &n);
	for (i=0;i<n;i++) {
	  printf("Step 9-%da. Low wavenumber limit of new forbidden region %d: ", i+1, i+1);
	  scanf("%lf", &xr1[nx+i]);
	  printf("Step 9-%db. High wavenumber limit of new firbidden region %d: ", i+1, i+1);
	  scanf("%lf", &xr2[nx+i]);
	}
	nx += n;
	printf("Step 10. Do you want to subtract an epoxy spectrum (1 for yes, 0 for no): ");
    scanf("%d", &Zflag);
	if (Zflag) {
	  printf("Step 10a. What is the filename of the backgrounded epoxy spectrum csv: ");
	  scanf("%s", epoxyfile);
	  printf("Step 10b. What wavenumber do you want to use to determine the expoy scaling (0 = default 3925.126): ");
	  scanf("%lf", &epoxypeakx); if (epoxypeakx == 0) epoxypeakx = 3925.126;
	}
	for (q=1,n=0;q<argc;q++) n+=arglist[q];
	printf("Step 11. How many input files not already given on command line (now there are %d; enter 0 for no more): ", n);
	scanf("%d", &zarg);
	for (i=0;i<zarg;i++) arglist[argc+i] = 1;
  }

  for (q=1,no=0;q<argc+zarg;q++) no += arglist[q];
  for (q=1;q<argc+zarg;q++) {
    if (arglist[q]) {
	  if (q < argc) strcpy(file1, argv[q]);
	  else {
	    printf("Next input filename: ");
		scanf("%s", file1);
	  }
      if ((fp1 = fopen(file1, "r")) == NULL) {
        printf("Could not open file %s\n", file1);
      } else {
	    printf("Reading file %s\n", file1);
        if (!o) {
		  n = 0;
          while (paulgetline(line, 100, fp1) != NULL) {
            n++;
          }
          rewind(fp1);
          x = (double *) malloc((unsigned) n*sizeof(double));
          spec = (double *) malloc((unsigned) n*sizeof(double));
          back = (double *) malloc((unsigned) n*sizeof(double));
		  iback = (double **) dmatrix(0,n,0,no);
		  if (oflag || oiflag) {
		    if (oflag) oback = (double **) dmatrix(0,n,0,no);
		    *allfiles = '\0';
		  }
		}
		if (oflag || oiflag) {
		  strcat(allfiles, ","); strcat(allfiles, file1);
		}
        for (i=0;i<n;i++) {
          paulgetline(line, 100, fp1);
	      sscanf(line, "%lf,%lf", &x[i], &spec[i]);
		  iback[i][o] = spec[i];
		  if (i==0) {
		    if (!o) x00 = x[0];
			else if (x[0] != x00) {
			  printf("WARNING: The first wavenumber in spectrum %s does not match the first spectrum...\n",file1);
			  printf("Only use spectra together that are sampled at exactly the same wavenumbers.\n");
			  printf("The program will continue now, but the results are not guaranteed.\n");
			}
		  }
        }
        fclose(fp1);
        strcpy(file2, file1);
		file2[strlen(file2)-3] = '\0';
		if (bgflag) strcat(file2, "bgs.csv");
		else strcat(file2, "bg.csv");
        if ((fp2 = fopen(file2, "w")) == NULL) {
          printf("Could not write file %s\n", file2);
        } else {
		  printf("Writing file %s\n", file2);
          n1 = 0;
          while (x[n1] < x1) n1++;
          n2 = 0;
          while (x[n2] < x2) n2++;
  
          /* First, 10 smoothing cycles */
          for (z=0;z<100;z++) {
            for (i=2; i<n-3; i++) {
	          spec[i] = (-spec[i-2] + 4.0*spec[i-1] + 4.0*spec[i+1] - spec[i+2])/6.0;
	        }
          }
		  
		  /*pre-zapping stage*/
		  for (j=0;j<npz; j++) {
	        pz1 = 0; while (x[pz1] < pzap1[j]) pz1++;
		    pz2 = 0; while (x[pz2] < pzap2[j]) pz2++;
		    for (i=pz1+1;i<pz2;i++) spec[i] = (spec[pz1]+(spec[pz2]-spec[pz1])/(x[pz2]-x[pz1])*(x[i]-x[pz1]));
          }
		  
          /* Now, 15-sample moving average */
          for (i=7;i<n-7;i++)
            for (j=-7,back[i]=0.0;j<=7;j++) back[i] += spec[i+j];
          for (i=7;i<n-7;i++) spec[i] = back[i]/15.0;
  
          for (i=0;i<n1;i++) back[i] = spec[i];
          for (i=n2;i<n;i++) back[i] = spec[i];
 
          /* Now locate anchor points */
          a[1] = n1; a[2] = n2;
		  for (z=3;z<3+nf;z++) {
		    a[z] = 0; while (x[a[z]] < f[z-3]) a[z]++;
			for (j=1;j<z;j++) {
	          if (a[z] < a[j]) {
	            swap = a[z];
		        for (i=z;i>j;i--) a[i] = a[i-1];
		        a[j] = swap;
		        break;
	          }
	        }
          }
          for (z=3+nf;z<=nanchors/*-ng*/;z++) {
            for (j=1;j<z-1;j++)
              for (i=a[j]; i<a[j+1]; i++)
                back[i] = spec[a[j]] + (spec[a[j+1]]-spec[a[j]])/(x[a[j+1]]-x[a[j]])*(x[i]-x[a[j]]);
            for (i=n1+1,max=0.0; i<n2; i++) {
			  for (xed=0,j=0;j<nx;j++) {
			    if (x[i] >= xr1[j] && x[i] <= xr2[j]) {
				  xed = 1;
				  break;
				}
			  }
              if (!xed && (back[i] - spec[i]) > max) {
	            max = back[i] - spec[i];
	        	a[z] = i;
	          }
	        }
	        if (max == 0.0) {
	          break;
	        }
	        for (j=1;j<z;j++) {
	          if (a[z] < a[j]) {
	            swap = a[z];
		        for (i=z;i>j;i--) a[i] = a[i-1];
		        a[j] = swap;
		        break;
	          }
	        }
          }
		  z--;
			
          /* Fill arrays */
          for (i=1;i<=z;i++) {
            xray[i] = x[a[i]];
	        y[i] = spec[a[i]];
          }
		  
		  if (!pflag) {
            /* Spline and splint */
            spline(xray, y, z, 1e32, 1e32, yy);
			for (j=1;j<z;j++) {
			  for (i=0,gflag=0;i<ng;i++) if (xray[j] < g[i] && g[i] < xray[j+1]) gflag = 1;
              for (i=a[j];i<a[j+1];i++) {
				if (gflag) back[i] = (y[j+1]*(x[i]-xray[j])+y[j]*(xray[j+1]-x[i]))/(xray[j+1]-xray[j]);
				else splint(xray, y, yy, z, x[i], &back[i]);
			  }
			}
		  } else {
		    double *w, *coeffs, **v, **a;
			int *yes;
			w = (double *) malloc((unsigned) (pdegree+1)*sizeof(double));
			coeffs = (double *) malloc((unsigned) (pdegree+1)*sizeof(double));
			v = (double **) dmatrix(1,pdegree+1,1,pdegree+1);
			a = (double **) dmatrix(1,z,1,pdegree+1);
			yes = (int *) malloc((unsigned) (z+1)*sizeof(int));
            for (i=1; i<=z; i++) yes[i] = 1;

			for (j=0;j<=pdegree;j++)
			  for (i=1;i<=z;i++) a[i][j+1] = pow(xray[i]/1000.0, (double) j);
			svdcmp(a, z, pdegree+1, w, v);
			for (j=1;j<=pdegree+1;j++) if (w[j] < 1.0e-08) w[j] = 0;
			svbksb(a, w, v, z, pdegree+1, y, yes, coeffs);

		    for (i=n1+1;i<n2;i++)
			  for (back[i]=0.0,j=0;j<=pdegree;j++) back[i] += coeffs[j+1]*pow(x[i]/1000.0, (double) j);
			
			free_dmatrix(a,1,z,1,pdegree+1);
			free_dmatrix(v,1,pdegree+1,1,pdegree+1);
			free(yes);
			free(coeffs);
			free(w);
		  }
		  
		  /*zapping stage*/
		  for (j=0;j<nz;j++) {
	        z1 = 0; while (x[z1] < zap1[j]) z1++;
		    z2 = 0; while (x[z2] < zap2[j]) z2++;
		    for (i=z1+1;i<z2;i++) back[i] += spec[i] - (spec[z1]+(spec[z2]-spec[z1])/(x[z2]-x[z1])*(x[i]-x[z1]));
          }
		  
		  /* epoxy subtraction stage */
		  if (Zflag) {
		    double clomp;
		    if (!o) {
			  if ((fp1 = fopen(epoxyfile, "r")) == NULL) {
			    printf("Could not open epoxy file, skipping epoxy subtraction.\n");
				Zflag = 0;
			  } else {
			    ne = 0;
                while (paulgetline(line, 100, fp1) != NULL) ne++;
                rewind(fp1);
                epoxyx = (double *) malloc((unsigned) ne*sizeof(double));
                epoxyy = (double *) malloc((unsigned) ne*sizeof(double));
			    for (i=0;i<ne;i++) {
			      paulgetline(line, 100, fp1);
				  sscanf(line, "%lf,%lf", &epoxyx[i], &epoxyy[i]);
			    }
			    for (i=0,clomp=10000;i<ne;i++) {
				  if (fabs(epoxyx[i] - epoxypeakx) < clomp) {
				    epoxypeakn = i; clomp = fabs(epoxyx[i] - epoxypeakx);
				  } else break;
				}
			  }
			}
			for (i=0,clomp=10000; i<n;i++) {
			  if (fabs(x[i] - epoxypeakx) < clomp) {
			   clomp = fabs(x[i] - epoxypeakx);
			  } else {
			    i-=1;
			    epoxyscale = (iback[i][o]-back[i])/epoxyy[epoxypeakn];
				break;
			  }
			}
			for (i=n1,j=0;i<n2;i++) {
			  while (epoxyx[j] < x[i]) j++;
			  back[i] += epoxyscale*epoxyy[j];
			}
		  }
			  
          if (bgflag) {
		    for (i=0;i<n;i++) {
			  fprintf(fp2, "%g,%g\n", x[i], iback[i][o]-back[i]);
			  if (oflag) oback[i][o] = iback[i][o] - back[i];
			}
		  } else {
		    for (i=0;i<n;i++) {
			  fprintf(fp2, "%g,%g\n", x[i], back[i]);
			  if (oflag) oback[i][o] = back[i];
			}
		  }
          fclose(fp2);
		  o++;
        }
		
		if (oiflag) {
		  if ((fp2 = fopen(ifile, "w")) == NULL) {
		    printf("Could not open omnibus input file %s, oops.\n", ifile);
		  } else {
		    fprintf(fp2, "%s\n", allfiles); 
		    for (i=0;i<n;i++) {
			  fprintf(fp2, "%g", x[i]);
			  for (j=0;j<no;j++) fprintf(fp2, ",%g", iback[i][j]);
			  fprintf(fp2, "\n");
			}
			fclose(fp2);
		  }
		}

		if (oflag) {
		  if ((fp2 = fopen(ofile, "w")) == NULL) {
		    printf("Could not open omnibus file %s, oops.\n", ofile);
		  } else {
		    fprintf(fp2, "%s\n", allfiles); 
		    for (i=0;i<n;i++) {
			  fprintf(fp2, "%g", x[i]);
			  for (j=0;j<no;j++) fprintf(fp2, ",%g", oback[i][j]);
			  fprintf(fp2, "\n");
			}
			fclose(fp2);
		  }
		}
	  }
	}
  }
}

void spline(double x[], double y[], int n, double yp1, double ypn, double y2[]) {
  int i, k;
  double p, qn, sig, un, *u;

  u = (double *) malloc((unsigned) n*sizeof(double));
  if (yp1 > 0.99e30)
    y2[1]=u[1]=0.0;
  else {
    y2[1] = -0.5;
    u[1] = (3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
  }
  for (i=2;i<=n-1;i++) {
    sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
    p=sig*y2[i-1]+2.0;
    y2[i]=(sig-1.0)/p;
    u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
    u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
  }
  if (ypn > 0.99e30)
    qn=un=0.0;
  else {
    qn=0.5;
    un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
  }
  y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
  for (k=n-1;k>=1;k--)
    y2[k]=y2[k]*y2[k+1]+u[k];
  free(u);
}

void splint(double xa[], double ya[], double y2a[], int n, double x, double *y) {
  static int klo=1, khi=1000;
  int k;
  double h, b, a;

  if (khi == 1000) khi = n;
  /* Allow increasing xa or decreasing xa */
  if (khi-klo != 1 || (xa[khi] < x && xa[klo] < x) || (xa[khi] > x && xa[klo] > x)) {
    klo=1;
    khi=n;
    while(khi-klo>1) {
      k=(khi+klo) >> 1;
      if ((xa[k] > x && xa[klo] <= x) || (xa[k] < x && xa[klo] >= x)) khi=k;
      else klo=k;
    }
  }
  h=xa[khi]-xa[klo];
  if (h==0.0) {
    printf("Splint failure...xa values must be distinct");
    return;
  }
  a = (xa[khi]-x)/h;
  b = (x-xa[klo])/h;
  *y = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}

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;
  }
}