#include <conio.h>
#include <math.h>
#include <graphics.h>

void   execute(void);
double aPoisson(unsigned n, float p, unsigned d);
long   factorial(unsigned x);

void main(void)
{
	char a;
	unsigned n,d;
	float p;

	clrscr();
	printf("This program plots operating characteristic (OC) curves \nof sampling plans using the poisson distribution.\n\n\tProgrammer: Ahmad Mostafa Gneady\n\n\n\n\n\n\n\n\n\nPress any key...");
	getch();

	do{
		execute();
		printf("\n\nAnother run?(y/n)");
		a=getche();
	}while(a=='y'||a=='Y');

	return;
}

void execute(void)
{
	char info[100];
	unsigned c,n;
	int drv=VGA,mode=VGAHI;
	float Pa,lq;

	clrscr();

	printf("Enter the sample size (n):");
	scanf("%u",&n);

	printf("\nEnter acceptance number (c)...");
	do{
		printf("\nc cannot be a large value (should be much smaller than %u):",n);
		scanf("%u",&c);
	}while(c>=n||c>12);
	sprintf(info,"sample size:%u, acceptance number:%u",n,c);

	/**************** START PLOTTING *******************/
	clrscr();
	initgraph(&drv,&mode,"");

	/******** AXES **********/
	setcolor(5);
	line(160,100,160,360);
	line(160,360,500,360);
	outtextxy(220,385,"lot quality (% defective)");
	setcolor(15);
	outtextxy(140,15,info);
	outtextxy(470,364,"100");
	outtextxy(156,364,"0");
	outtextxy(150,355,"0");
	outtextxy(150,115,"1");
	setcolor(9);
	outtextxy(30,420,"If you are running this program under Windows, you can");
	outtextxy(30,430,"press <Print Screen> to copy this screen to the Windows Clip Board.");
	outtextxy(30,450,"Press any key to continue...");
	settextstyle(0,1,1);
	setcolor(5);
	outtextxy(130,130,"Probability of acceptance");

	moveto(160,120);
	setcolor(15);
	for(lq=0.01;lq<1.00;lq+=.01)
	{
		Pa = aPoisson(n,lq,c);
		lineto(640*(lq+.5)/2,480-480*(Pa+.5)/2);
	}
	lineto(640*1.5/2,480-480*.5/2);
	getch();
	closegraph();

	return;
}

double aPoisson(unsigned n, float p, unsigned d)
{
	int i;
	double s=0.0;
	for(i=0;i<=d;i++)
		s+=pow(n*p,i)*exp(-1.0*n*p)/factorial(i);
	return s;
}

long factorial(unsigned x)
{
	int i;
	long y=1;

	if(!x) return 1;
	for(i=1;i<=x;i++)
		y*=i;
	return y;
}