#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>

/*Constants*/

#define	MAXSAMPLESIZE 25
#define	MAXOBS 50
#define FALSE 0
#define	TRUE !FALSE
#define	OK TRUE


int	datatype;		//data type 1 or 2
int	samplesize; 		//number of observed values per sample, n
int	numofsamples;		//number of observations collected, m
char	*filename = "temp.sqc";
FILE	*fn;
double	sample[MAXOBS+1][MAXSAMPLESIZE+1];	//input data, stores either sample data or number of defects
double	usl,lsl;		//specification limits, optional

double	rhi, rlow, rrange;	// variables for graphing charts



/* Control Chart Parameters */

static double A2[] = { 1.880, 1.023, 0.729, 0.577, 0.483, 0.419, 0.373, 0.337, 0.308, 0.285, 0.266, 
0.249, 0.235, 0.223, 0.212, 0.203, 0.194, 0.187, 0.180, 0.173, 0.167, 0.162, 0.157, 0.153};
static double d2[] = {1.128, 1.693, 2.059, 2.326, 2.534, 2.704, 2.847, 2.970, 3.078, 3.173, 3.258, 
3.336, 3.407, 3.472, 3.532, 3.588, 3.640, 3.689, 3.735, 3.778, 3.819, 3.858, 3.895, 3.931};
static double D3[] = { 0.000, 0.000, 0.000, 0.000, 0.000, 0.076, 0.136, 0.184, 0.223, 0.256, 0.283, 
0.307, 0.328, 0.347, 0.363, 0.378, 0.391, 0.403, 0.415, 0.425, 0.434, 0.443, 0.451, 0.459};
static double D4[] = { 3.267, 2.575, 2.282, 2.115, 2.004, 1.924, 1.864, 1.816, 1.777, 1.744, 1.717, 
1.693, 1.672, 1.653, 1.637, 1.622, 1.608, 1.597, 1.585, 1.575, 1.566, 1.557, 1.548, 1.541};


/* Functions prototypes */

void graphit(int, int, int, double*, double, double, double);
static read_init(int k);
static int read_additional(int k);
static double mean(int r, int k);
static double range(int r, int k);
static process_capability(double, double);

/* Main program */

int main(int argc, char *argv[])             
{
  double	samplemean[MAXOBS+1];
  double	samplerange[MAXOBS+1];  
  double	sampleproportion[MAXOBS+1];
  double	grandmean=0.0;
  double	avgrange=0.0;
  double	avgproportion=0.0;		
  int		i,j;
  char*		datatypename[] = {"", "Variable(measure)","Attribute"};
  double	cl_x,lcl_x,ucl_x;  /* center line, lower control limit, upper control limit for X-chart*/
  double	cl_r,lcl_r,ucl_r;  /* center line, lower control limit, upper control limit for R-chart*/
  double	cl_p,lcl_p,ucl_p;  /* center line, lower control limit, upper control limit for P-chart*/
  int		totalnum;	   /* total number of observations:numofsamples+additional */
  static char help[] = "chart datatype samplesize numofsamples [usl lsl]\n\
   datatype = 1-2, where 1=Variable(measure), 2=Attribute\n\
   samplesize = number of observed values per sample, <=25\n\
   numofsamples = number of observations collected\n\
   usl = upper specification limit, only for datatype=1\n\
   lsl = lower specification limit, only for datatype=1\n";



  if (argc<4)  {printf(help); exit(0);}

  datatype=atoi(argv[1]);
  samplesize=atoi(argv[2]);
  numofsamples=atoi(argv[3]);
  if (argc==6) {usl=atof(argv[4]);lsl=atof(argv[5]);}

  printf("\nStatistical Quality Control For %s Data Type\n\n",datatypename[datatype]);
  printf("Input values:\n");
  printf("------------------- \n");
  printf("Data type:\t\t%i, %s\nSample size, n: \t%d \nNumber of samples, m: \t%d\n",datatype,datatypename[datatype],samplesize,numofsamples);
    

  if (datatype==1)
	{
	if (samplesize>25)
	  {printf("\nSample size should be <=25\n\n");return 0;}  
	
	read_init(samplesize);			/* reading from the file model based sample data*/
	
	printf("\n\nSample\t\tObserved values\t\tMean\tRange\n");
        printf("------\t-----------------------\t\t-----\t------\n");
	
	for (i=1;i<=numofsamples;++i)
	  {
	  samplemean[i]=mean(i,samplesize);
	  grandmean=grandmean+samplemean[i];	/* calculating average and range for each sample */
	  samplerange[i]=range(i,samplesize);
	  avgrange=avgrange+samplerange[i];
	  }
	grandmean=grandmean/numofsamples;	/* calculating grand mean X-bar */
	avgrange=avgrange/numofsamples;		/* calculating average range R-bar */

	for (i=1;i<=numofsamples;++i)		/* output on the screen */
	  {
	  printf("%d\t",i);
	  for (j=1;j<=samplesize;++j)
	    {
	    printf("%4.1lf ",sample[i][j]);
	    }
	  if (samplesize<4) printf("\t\t"); if (samplesize==4) printf("\t");  
	  printf("\t%4.1lf\t%4.1lf",samplemean[i],samplerange[i]);
	  printf("\n");
	  }
	printf("\nGrand mean:\t%5.2lf\nAverage range:\t%5.2lf\n\n",grandmean,avgrange);
	
						//*calculating CL,UCL, and LCL
	if (samplesize >=2)
		{
		cl_x = grandmean;
		ucl_x = grandmean + A2[samplesize-2]*avgrange;	
		lcl_x = grandmean - A2[samplesize-2]*avgrange;	
		if (lcl_x<0) {lcl_x = 0;}
		cl_r = avgrange;
		ucl_r = D4[samplesize-2]*avgrange;
		lcl_r = D3[samplesize-2]*avgrange;
		printf("\t\tX-chart\tR-chart\n");
		printf("\t\t-------\t-------\n");
		printf("UCL\t\t%6.2lf\t%6.2lf\n",ucl_x,ucl_r);
		printf("CL \t\t%6.2lf\t%6.2lf\n",cl_x,cl_r);
		printf("LCL\t\t%6.2lf\t%6.2lf\n",lcl_x,lcl_r);

		totalnum=read_additional(samplesize);   
		for (i=numofsamples+1;i<=totalnum;++i)
		  {
		  samplemean[i]=mean(i,samplesize);	/* calculating sample mean and range for */
		  samplerange[i]=range(i,samplesize);	/* additional sample data entered */
		  }

		graphit(1,numofsamples,totalnum,samplemean,lcl_x, cl_x, ucl_x);		/* graphing x-chart and */
		graphit(2,numofsamples,totalnum,samplerange,lcl_r, cl_r, ucl_r);	/* r-chart */

		/* process capability study if usl&lsl are given */

		if ((usl!=0)&&(usl>=lsl))
		  {
		  process_capability(grandmean,avgrange);
		  }
	  	} // end if (sample size>=2)
	else {printf("Sample size is too small to calculate characteristics for X-chart and R-chart\n\n");}
	} // end if (datatype==1)

  if (datatype == 2) 
	{
	read_init(1);			/* reading sample data from the file */
		
	for (i=1;i<=numofsamples;++i)
	  {
	  sampleproportion[i]=sample[i][1]/samplesize;
	  /* check if input data is correct */
	  if (sampleproportion[i]>1)
	    {
		printf("\nNumber of defective cannot be more than sample size.\nInput correct data.\n");
		return 0;
	    }
	  avgproportion=avgproportion+sampleproportion[i];
	  }
	avgproportion=avgproportion/numofsamples;

	printf("\n\nSample\t# of defects\tFraction defective\n");
        printf("------\t------------\t-----------------\n");

	/* check if sample size is large enough */
		
	if ((avgproportion*samplesize>=5)||(samplesize*(1-avgproportion))) 
	{
		for (i=1;i<=numofsamples;++i)
		  {
		  printf("%d\t",i);
		  printf("%4.0lf\t ",sample[i][1]);
		  printf("\t%4.2f",sampleproportion[i]);
		  printf("\n");
	  	}
		printf("\nAverage proportion:\t%6.4lf\n\n",avgproportion);
	
							// *calculating CL,UCL, and LCL *//

		cl_p = avgproportion;
		ucl_p = avgproportion + 3*sqrt(avgproportion*(1-avgproportion)/samplesize);	
		lcl_p = avgproportion - 3*sqrt(avgproportion*(1-avgproportion)/samplesize);	
		if (lcl_p<0) {lcl_p = 0;}
		printf("\t\tP-chart\n");		// * output *//
		printf("\t\t-------\n");
		printf("UCL\t\t%6.4lf\n",ucl_p);
		printf("CL \t\t%6.4lf\n",cl_p);
		printf("LCL\t\t%6.4lf\n",lcl_p);

		totalnum = read_additional(1);	/* reading additional sample data from file for display only */

		for (i=numofsamples+1;i<=totalnum;++i)
		  {
		  sampleproportion[i]=sample[i][1]/samplesize;
		  /* check if input data is correct */
		  if (sampleproportion[i]>1)
		    {
			printf("\nNumber of defective cannot be more than sample size.\nInput correct data.\n");
			return 0;
		    }
		  }
		graphit(3,numofsamples,totalnum,sampleproportion,lcl_p, cl_p, ucl_p);	/* graphing p-chart */
  	}
	  else	printf("\nSample size is not large enough. Increase it and repeat");
	
	}//end if (datatype==2)
  
//end of main
return 0;
}

/* Reading from the file to the array, model-based data */
read_init(k)
int	k;	// number of values per sample to read from the file (either sample size or 1)
{
	int i,j;
	double	q;
	int nscan;	// to detect end of file
	
	fn=fopen(filename,"r");
	for (i=1;i<=numofsamples;++i)
	{
	   for (j=1;j<=k;++j)
	   {
		nscan=scanf ("%lf",&q);
		if (nscan==EOF) 
		   {
		    printf("\nThere is not enough sample data\nbased on entered sample size or number of observations.\n");
		    exit(0);
		   }
		sample[i][j] = q;
	   }
        }
}

/* Reading from the file to the array, additional data used for display only */
int read_additional(k)
int	k;	// number of values per sample to read from the file (either sample size or 1)
{
	int i,j;
	double	q;
	int nscan;	// to detect end of file
	
	nscan=scanf ("%lf",&q);	
	i=numofsamples+1;
	while (nscan!=EOF)
	{
	   for (j=1;j<=k;++j)
	   {
		sample[i][j] = q;
		nscan=scanf ("%lf",&q);
				
	   }
	   ++i;
        }
	return(i-1);	// returns total number of entered observations 
}

/* Finding sample mean */
double mean(r,k)
int r;
int k;
{
	int j;
	double sum;

	sum=0.0;
	for (j=1;j<=k;++j)
	{
	  sum=sum+sample[r][j];
	}	
	sum=sum/k;
	return sum;
}


/* Finding sample range */
double range(r,k)
int r;
int k;
{
	int j;
	double diff;
	double max,min;

	diff=0.0;
	max=sample[r][1];
	min = max;
	for (j=2;j<=k;++j)
	{
	   if (sample[r][j]>max) max=sample[r][j];
	   if (sample[r][j]<min) min=sample[r][j];   
	}
	diff=max-min;
	return diff;
}

/* Analyzing Process Capability if specifications are given */
process_capability (double x, double r)
{
	double sigma;
	double cp, cpk, sigma_cap_index;
	double cpk1,cpk2;

	sigma = r/d2[samplesize-2];	/* process standard deviation */
	cp=(usl-lsl)/6/sigma;		/* process capability ratio */
	cpk1=(usl-x)/3/sigma;
	cpk2=(x-lsl)/3/sigma;
	cpk=cpk1;
	if (cpk2<cpk1) cpk=cpk2;
	sigma_cap_index = 1.5+3*cpk;
	printf("\n\nProcess Capability:\n-------------------\nGiven Spec Limits %7.4lf and %7.4lf:\n\n",usl,lsl);
	printf("Process standard deviation\t%5.2lf\nCp\t\t\t\t%7.4lf\nCpk=min{%7.4lf,%7.4lf}\t%7.4lf\n",sigma,cp,cpk1,cpk2,cpk);
	printf("Sigma Capability Index\t\t%5.2lf\n\n",sigma_cap_index);

}

static inline int graphmark(double value) { return (int) 1+15*2*(value-rlow)/rrange; }

void graphit(int ctype, int m, int m2, double *x, double lcl, double cl, double ucl)
/*
 *  ctype = chart type: 1=xbar, 2=R, 3=p
 *  m = number of statistical values
 *  m2 = number of statistical values, including additional observations
 *  x[] = vector of m statistics, with first value in x[0]
 *  lcl = lower control limit
 *  cl  = center line
 *  ucl = upper control limit
 */
	{
	register int i,j;
	int ilcl, icl, iucl, mdiff;
	double	res;
	static char *fhd     = "%4s %9s[i]%12s  %s\n";
	static char *fhd2    = "%-25s      %s\n";
	static char buf[]     =  "                                 ";
	static char heading[] =  "---------------------------------";
	static char type[][20] = { "x-bar", "R", "p"};
	char isave;
	char *criteria;       // Out of control reason or blank
	int ook;    // logical, out of control?
	static char *blanks = " ";
	int run;    // count of # times on one side of CL
	int side;   // which side is on
	int abovemean;  // logical for "is above mean?"

	// normalize data first and calculate heading strings here
	
	printf("\n\t\t\t*** %s CONTROL CHART  ***\n",&type[ctype-1][0]);
	printf("\nLCL = %lf   CL = %lf   UCL = %lf\n\n",lcl,cl,ucl);

	for (i=1,rhi=ucl,rlow=lcl; i <= m2; ++i) 
			{
			if (x[i] > rhi) rhi = x[i];
			if (x[i] < rlow) rlow = x[i];
			}
	rrange = rhi-rlow;
	ilcl = graphmark(lcl);
	icl  = graphmark(cl);
	iucl = graphmark(ucl);

	heading[icl] 	= 'C';
	heading[ilcl]	= 'L';
	heading[iucl]	= 'U';
	buf[icl] = buf[ilcl] = buf[iucl] = '|';
	mdiff = (m != m2);

	printf(fhd,"i.",&type[ctype-1][0],"   ",heading);

	abovemean = (x[1]>cl);
	run = 0;
	for (i=1; i<=m2; ++i)
		{
		criteria = blanks;
		ook = FALSE;
		side = (x[i]>cl);
		if(side == abovemean) {++run;}
		else {abovemean = side; run=1;}
		if (run >= 6) {criteria = " Run 6+"; ook = TRUE;}
		if (x[i] > ucl) {criteria = " > UCL"; ook = TRUE;}
		if (x[i] < lcl) {criteria = " < LCL"; ook = TRUE;}

		printf("%3d. %12.2lf%-12s", i,x[i],criteria);
		j		= graphmark(x[i]);
		isave   = buf[j];
		buf[j]	= (ook)? '#' : '*';
		printf("  %s\n",buf);
		buf[j]	= isave;
		if (mdiff && (i == m)) printf(fhd2,"Additional observations:",heading);
		}
	printf(fhd,"",&type[ctype-1][0],"",heading);
	buf[icl] = buf[ilcl] = buf[iucl] = ' ';
	heading[icl] =  heading[ilcl] = heading[iucl] = '-';
	}