#include <stdio.h>
#include <graphics.h>
#include <conio.h>
#include <math.h>
#include <stdlib.h>
#include <dos.h>

 int ip,n,i,j,k,ns,nn,j1,i1,nunk,jj, nx, ny, ni;
 double ch, energy, d, offset,dele,delch,dele2;
 long int indx[10],nsi[10];
 double col[10],xcal[10],fit[10],delfit[10],ycal[10];
 double y[10],x[10],erry[10];
 int running;
 int graphdriver, graphmode, error;
 int maxx,maxy,yp;
 char cho;
 double xunk[10], errunk[10];
 float sumx,sumy,sumxy,sumx2,sumx_sq,slope,yfit,intercept,oset;
 float sumdif,xscale,yscale,xmin,ymin,xmax,ymax,xi;
 float delyy, delm, delyi, delxi, dett;

void linfit();
void graphit();
void screenin();
void readfl();


void main()
{

printf(" Welcome to the linear fitting program \n\n");

running=1;
printf(" press s to input data from the screen\n");
printf(" press r to read data from a file\n");
do
 {
  if (kbhit)
   {
     cho=getch();
     switch(cho)
     {
      case 's' : screenin();
		 running=2;
		 break;
      case 'r' : readfl();
		 running=2;
		 break;
      default :  printf("Try again, s or r\n");
		 break;
     }
 }
}while (running==1);
running=1;
printf("\n\n");
printf(" press l to do a least-square linear fit \n");
printf(" press q to exit the program\n");

do
 {
  if (kbhit)
   {
     cho=getch();
     switch(cho)
     {
      case 'q' : closegraph();
		 running=2;
		 break;

      case 'l' : closegraph();
		 linfit();
		 graphit();
		 break;

      default :  closegraph();
		 printf("Try again, q or l\n");
		 break;
     }
  }
}  while (running==1);
}

void linfit()
{
printf("\n\n");
printf(" number       x-value       y-value  \n");
for (i = 1; i<= nn; i++)
printf("   %d         %8.3f         %8.3f   \n", i, x[i],y[i]);
printf(" How many points from the above list do you \n");
printf(" want to use for your fit? \n");
scanf("%d",&ns);
printf(" Which ones?  Enter the numbers in order\n");
for (i=1; i<=nn; i++)
      nsi[i]=0;
for (i=1; i<=ns; i++)
      scanf("%d",&nsi[i]);

for (i=1; i<=ns; i++)
 {
  xcal[i]=x[nsi[i]];
  ycal[i]=y[nsi[i]];
 }
sumx=0;
sumy=0;
sumx2=0;
sumxy=0;
for (i=1; i<=ns; i++)
 {
  sumx=sumx+xcal[i];
  sumy=sumy+ycal[i];
  sumx2=sumx2+xcal[i]*xcal[i];
  sumxy=sumxy+xcal[i]*ycal[i];
 }
sumx_sq=sumx*sumx;
dett=ns*sumx2-sumx_sq;
slope=(ns*sumxy-sumx*sumy)/(ns*sumx2-sumx_sq);
intercept=(sumy*sumx2-sumxy*sumx)/(ns*sumx2-sumx_sq);
oset=-intercept/slope;
sumdif=0.0;
for (i=1; i<=ns; i++)
 {
  sumdif=sumdif+pow(slope*xcal[i]+intercept-ycal[i],2);
 }
delyy=sqrt(sumdif/ns);
delyi=delyy*sqrt(sumx2/dett);
delm=delyy*sqrt(ns/dett);
delxi=oset*(delm/slope+delyi/intercept);
if (delxi<0) delxi=-delxi;
xmax=x[1];
for (i=2; i<=nn; i++)
 if (x[i]>xmax) xmax=x[i];
ymax=y[1];
for (i=2; i<=nn; i++)
 if (y[i]>ymax) ymax=y[i];
xmin=x[1];
for (i=2; i<=nn; i++)
 if (x[i]<xmin) xmin=x[i];
ymin=y[1];
for (i=2; i<=nn; i++)
 if (y[i]<ymin) ymin=y[i];
xmin=0.0;
ymin=0.0;
}

void graphit()
{
char buffer[80];

graphdriver=DETECT;
initgraph(&graphdriver,&graphmode," ");
maxx=getmaxx();
maxy=getmaxy();
yp=340;
line(5,0,5,300);
line(5,300,maxx,300);
xscale=.95*(maxx-5.0)/(xmax-xmin);
yscale=295.0/(ymax-ymin);
sprintf(buffer,"Slope = %10.4f   y-int = %10.4f   x-int = %10.4f ",slope,intercept,oset);
outtextxy(20,yp-20,buffer);
sprintf(buffer,"Del m = %10.4f  del yi = %10.4f  del xi = %10.4f ",delm,delyi,delxi);
outtextxy(20,yp-10,buffer);
ip=1;
yp=yp+10;
outtextxy(20,yp,"Standards:   x-value     y-value        Difference ");
for (jj=1; jj<=nn; jj++)
{
  yfit=slope*x[jj]+intercept;
  sprintf(buffer,"%15.2f    %15.2f    %15.2f",x[jj],yfit,y[jj]-yfit);
  if (jj==nsi[ip])
   {
    setcolor(13);
    outtextxy(70,floor(yp+jj*10),buffer);
    nx=floor(5+(x[jj]-xmin)*xscale);
    ny=floor(300-(y[jj]-ymin)*yscale);
    circle(nx,ny,3);
    ip=ip+1;
   }
    else
     {
      setcolor(10);
      outtextxy(70,floor(yp+jj*10),buffer);
      nx=floor(5+(x[jj]-xmin)*xscale);
      ny=floor(300-(y[jj]-ymin)*yscale);
      circle(nx,ny,3);
     }
 }
yp=yp+jj*10+20;

setcolor(14);
outtextxy(20,yp,"Unknowns:   x-value        y-fit");
settextstyle(3,0,1);

for (i=1; i<=nunk; i++)
 {
  yfit=slope*xunk[i]+intercept;
  sprintf(buffer,"%15.2f     %15.2f",xunk[i],yfit);
  outtextxy(60,yp+i*15,buffer);
  nx=floor(5+(xunk[i]-xmin)*xscale);
  ny=floor(300-(yfit-ymin)*yscale);
  circle(nx,ny,3);
 }
 for (i=0; i<=maxx; i++)
  {
   xi=i/xscale+xmin;
   ny=floor((slope*xi+intercept - ymin)*yscale);
   putpixel(floor(5+i),floor(300-ny),13);
  }
 setcolor(15);
 settextstyle(0,0,0);
 outtextxy(10,maxy-10," l: for least-sq fit;           q: quit");

}

void screenin()
{
 char filename[80];
 FILE *stream;

 printf(" input the number of data points (limit=10)\n");
 scanf("%d",&nn);
for (i=1; i<=nn; i++)
 {
  printf(" input the x-values and y-values for data point %d\n",i);
  scanf("%lf  %lf",&x[i],&y[i]);
  printf("\n");
 }
printf("input the number of interpolation values (limit=10)\n");
scanf("%d",&nunk);
for (i=1; i<=nunk; i++)
 {
  printf(" input the x-value for the unknown %d \n",i);
  scanf("%lf",&xunk[i]);
 }
printf("filename to save file: \n");
scanf("%s", &filename);
stream=fopen(filename,"w");
fprintf(stream,"%d %d \n",nn,nunk);
for (i=1; i<=nn; i++)
  fprintf(stream,"%lf    %lf\n",x[i],y[i]);
for (i=1; i<=nunk; i++)
  fprintf(stream,"%lf \n",xunk[i]);
fclose(stream);
return;
}

void readfl()
{
 char filename[80];
 FILE *stream;

 printf("\n");
 printf("file to be read in:\n");
 scanf("%s",&filename);
 stream=fopen(filename,"r");
 fscanf(stream,"%d %d",&nn,&nunk);
for (i=1; i<=nn; i++)
  fscanf(stream,"%lf %lf",&x[i],&y[i]);
for (i=1; i<=nunk; i++)
  fscanf(stream,"%lf",&xunk[i]);
fclose(stream);
printf("\n");
printf("Calibration energys, channel number \n");
for (i=1; i<=nn; i++)
 printf("%lf    %lf \n",x[i],y[i]);
printf("\n Unknown channel numbers \n");
for (i=1; i<=nunk; i++)
 printf("%lf \n",xunk[i]);
return;
}