#include <utility.h>
#include <ansi_c.h>
#include <stdlib.h>
#include <cvirte.h>		
#include <userint.h>
#include "assign4_cvi.h"

#define ON  1
#define OFF 0
#define YES 1
#define NO  0

int runtrancalc(int num_cell, double delay_time, int num_step, double dt);
int runsscalc(int num_cell, double delay_time, int num_iter);
void plotdata(void);

static int panelHandle, y1Handle;  

float trace[2]; //strip chart input



int main (int argc, char *argv[])
{
	if (InitCVIRTE (0, argv, 0) == 0)
		return -1;	/* out of memory */
	if ((panelHandle = LoadPanel (0, "assign4_cvi.uir", PANEL)) < 0)
		return -1;
	DisplayPanel (panelHandle);
	RunUserInterface ();
	DiscardPanel (panelHandle);
	return 0;
}

int CVICALLBACK QuitCallback (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
	switch (event)
		{
		case EVENT_COMMIT:
			QuitUserInterface(0);
			break;
		}
	return 0;
}

int CVICALLBACK RunTransientCallback (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
	int size, num_step, num_cell;
	double delay_time, dt;
	switch (event)
		{
		case EVENT_COMMIT:
			GetCtrlVal (panelHandle, PANEL_DT, &dt);
			GetCtrlVal (panelHandle, PANEL_DELAY, &delay_time);
			GetCtrlVal (panelHandle, PANEL_NUMSTEPS, &num_step);
			GetCtrlVal (panelHandle, PANEL_NUMCELL, &num_cell);
			size = runtrancalc(num_cell, delay_time, num_step, dt);
			//plotdata(size, flux);
			break;
		}
	return 0;
}

int CVICALLBACK RunSSCallback (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
	int size, num_iter, num_cell;
	double delay_time;
	switch (event)
		{
		case EVENT_COMMIT:
			GetCtrlVal (panelHandle, PANEL_DELAY, &delay_time);
			GetCtrlVal (panelHandle, PANEL_NUMITER, &num_iter);
			GetCtrlVal (panelHandle, PANEL_NUMCELL, &num_cell);
			size = runsscalc(num_cell, delay_time, num_iter);
			break;
		}
	return 0;
}


int CVICALLBACK ResetCallback (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
	switch (event)
		{
		case EVENT_COMMIT:
			DeleteGraphPlot (panelHandle, PANEL_GRAPH, -1, VAL_IMMEDIATE_DRAW);  
			ClearStripChart (panelHandle, PANEL_STRIPCHART);
			break;
		}
	return 0;
}

int runtrancalc(int num_cell, double delay_time, int num_step, double dt)

// EP4D3 Assignment 4, question 3  solution   2003.11.14
// Case: Slab reactor, transient
// Method: Semi-implicit step

{
 	int i, j, forward, n;
 	float t=0.0, tfin, v=??;
 	float x[1000+2],flux[1000+2], flux_old[1000+2];  //make room for real cells + phantoms at the edges
                  //size = num_cell+2
 	float s[1000+2];

 	float h;  //cell size
 	float D=??  // diffusion coefficient
 	float ea =??   //absorption cross section
 	float nu =??  //nu in the fission term
 	float ef =?? // sigma f
 	float k =?? // k fudge factor
 	float max_flux, s_old, s_new;  //tmp value holders
 	float a = ??, b = ??, flux_m, flux_m_old; //controller parameters

	char output_filename[13]="tran.dat";
 	FILE *fp;

 	h =??;  //cell size

 	if (( fp = fopen(output_filename,"w")) == NULL)
 	{
   		printf("Could not open %s",output_filename);
   		return(0);
 	}
 	fprintf(fp, "Transient results\n");

 	tfin=dt*num_step;

 	x[0] =??;
	s_old = 0.0;
 	for (i=1;i<=num_cell;i++)
 	{
   		flux[i]=flux_old[i]=1.0;   // initial guess at flux
   		x[i] = x[i-1]+h; // calc mesh distance
 	}

 
 	flux[0]=flux_old[0]=??;
 	flux[num_cell+1]=flux_old[num_cell+1]=??;      // flux is ?? at edges
 	x[num_cell+1] = x[num_cell]+h;
 	
 	forward = 1;
 	while (t<=tfin)
 	{
   		t=t+dt;
   		//printf("\n%.8f",t);
   		fprintf(fp,"\n%.8f",t);
   		/* this is the equation derived from the difference formula */
   		//sweep over all space, left to right on one pass, then right to left on the next pass
   		for (n=1;n<=num_cell;n++)
   		{
     		if (forward ==1) i=n; //normal forward sweep
     		else i=num_cell+1-n; //reverse sweep
     		s[i] = ??; // set the source as ??
 	 		flux_old[i] = flux[i];
     		flux[i]=??
     		//printf("\t%.4f",flux[i]);
     		fprintf(fp,"\t%.4f",flux[i]);
   		}
   		//change the sweep direction parameters to alternate direction on next sweep
   		if (forward == 1) forward = 0;
   		else forward = 1;

		//update k
		
		flux_m = flux[(num_cell+1)/2];
		flux_m_old = flux_old[(num_cell+1)/2];
		k = k + a * (flux_m - 1.0) + b * ( flux_m - flux_m_old)/dt;
 
		SetCtrlVal (panelHandle, PANEL_K, k);

   		y1Handle = PlotXY (panelHandle, PANEL_GRAPH, x, flux, num_cell+2,
					   VAL_FLOAT, VAL_FLOAT, VAL_THIN_LINE,
					   VAL_SOLID_DIAMOND, VAL_SOLID, 1, VAL_RED);
		trace[0]=k;
		trace[1]=flux_m;
		PlotStripChart (panelHandle, PANEL_STRIPCHART, trace, 2, 0, 0, VAL_FLOAT);
		Delay (delay_time);
  	}

 fclose( fp );
return (num_cell);
}
//---------------------------------------------------------------------------

int runsscalc(int num_cell, double delay_time, int num_iter)

// EP4D3 Assignment 4, question 3  solution   2003.12.14
// Case: Slab reactor, steady state
// Method: Semi-implicit step, combined inner and outer loops

{
	int i, j, forward, n;
	float x[1000+2],flux[1000+2], flux_old[1000+2];  //make room for real cells + phantoms at the edges
    			               //size used = num_cell+2
 	float s[1000+2];	//neutron source 
 	float h;  //cell size
 	float D=??  // diffusion coefficient
 	float ea =??   //absorption cross section
 	float nu =??  //nu in the fission term
 	float ef = // sigma f
 	float k =1.0; // k fudge factor
 	float max_flux, s_old, s_new;  //tmp value holders
 
 	char output_filename[13]="steady.dat";
 	FILE *fp;

 	h =?  //cell size
 	
 	if (( fp = fopen(output_filename,"w")) == NULL)
 	{
   		printf("Could not open %s",output_filename);
   		return(0);
 	}
 	fprintf(fp, "Steady State results\n");

 	x[0] = -60.0;
	s_old = 0.0;
 	for (i=1;i<=num_cell;i++)
 	{
   		flux[i]=flux_old[i]=1.0;   // initial guess at flux
   		x[i] = x[i-1]+h; // calc mesh distance
 	}

 	flux[0]=flux_old[0]=??
 	flux[num_cell+1]=flux_old[num_cell+1]=??      // flux is ?? at edges
 	x[num_cell+1] = x[num_cell]+h;
 	
	//initialize sweep control parameters for a left to right sweep
	forward = 1;
 	//iterate to convergence or till # of iterations reached
 	for (j=0;j<num_iter;j++)
 	{
   		fprintf(fp,"\n%d ",j);
   		//sweep over all space, left to right on one pass, then right to left on the next pass
   		for (n=1;n<=num_cell;n++)
   		{
     		if (forward ==1) i=n; //normal forward sweep
     		else i=num_cell+1-n; //reverse sweep
     		s[i] = ??
 	 		flux_old[i] = flux[i];
     		/* this is the equation derived from the difference formula */
	   		flux[i]=??
     		//printf("\t%.4f",flux[i]);
     		fprintf(fp,"\t%.4f",flux[i]);
   		} //end of sweep over all space
   		//change the sweep direction parameters to alternate direction on next sweep
   		if (forward == 1) forward = 0;
   		else forward = 1;
   		
		//update k
		s_new = 0.0;
		s_old = 0.0;
		//calc the source terms
		for (i=1; i<=num_cell; i++)
		{
			s_new = s_new + ??
			s_old = s_old + ??
		}
		k = ??
 
		SetCtrlVal (panelHandle, PANEL_K, k);


		//normalize the flux and save flux vector
		
		max_flux = flux[(num_cell+1)/2];
    	for (i=1;i<=num_cell;i++)
   		{
   			if (max_flux > 0.0) flux[i] = flux[i]/max_flux;
		}
		
   		y1Handle = PlotXY (panelHandle, PANEL_GRAPH, x, flux, num_cell+2,
					   VAL_FLOAT, VAL_FLOAT, VAL_THIN_LINE,
					   VAL_SOLID_DIAMOND, VAL_SOLID, 1, VAL_RED);
		trace[0]=k;
		trace[1]=max_flux;
		PlotStripChart (panelHandle, PANEL_STRIPCHART, trace, 2, 0, 0, VAL_FLOAT);
		Delay (delay_time);
		
  	}  // end of outer (source) iteration loop
  

 	fclose( fp );
	return (num_cell);
}
//---------------------------------------------------------------------------