#include <X11/Xlib.h>
#include <assert.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <sys/io.h>
#include <unistd.h>
#include <asm/msr.h>
#include <fcntl.h>
#include <errno.h>
#include <string.h>
#include <sys/mman.h>
#include <sys/ioctl.h>
//  #include <asm/page.h>
#include <sys/time.h>
#include <time.h>
#include <termios.h>
#include <sys/signal.h>
#include <sys/types.h>

#define BAUDRATE B115200
#define MODEMDEVICE "/dev/ttyUSB0"
#define _POSIX_SOURCE 1 /* POSIX compliant source */
#define FALSE 0
#define TRUE 1

volatile int STOP=FALSE;

int wait_flag=TRUE;                    /* TRUE while no signal received */


int xi,y;
int ppc_fix;
int ii,i, j, xn, yyn, icolor;
int xx1, xx2, yy1, yy2;
//  Global
int running, run, runc, graphdriver, graphmode, graphing;
int nx,ny, ncheck;
long int count;
float con1, con2, sumtime, v1c, v2c, xa, denom, circum;
float a, sumxx, b, sumxy, sumx, sumy, del, delx, dely;
float maxx, maxy, maxxx, maxxy, scalet, scalev;
char cho, filename[20], buffer[80], dum[80];
double on[100], off[100], fac, cpu, fr, frc,facpic,tottime;
float x, time1[100], vel[100], rvel[100], toffset;
int ntot, ncal, s1, s2;
float ss1, inter, slop;
unsigned char data[256];
unsigned long tl[100], th[100], ontl[100], onth[100], ptest;
unsigned long dum1, dum2, dum3, dum4;
FILE *inf;
int fd,c, res, icount;
struct termios oldtio,newtio;
struct sigaction saio;           /* definition of signal action */
unsigned char buf[255];

void menu();
void calbrat();
void setpoints();
void takedat();  //244
void datview();  //157
void velplot();  //167
void testgate();  //98
void signal_handler_IO (int status);   /* definition of signal handler */

int main()
{
  printf("\n\n\n   CONSTANT ACCELERATION PROGRAM \n\n\n");
  printf("Software by Rich Aguayo (Cal Poly Pomona Physics graduate 91)\n");
  printf("        Ray Chen (Cal Poly Pomona Physics graduate 04)\n");
  printf("   and Andrew Avila (Biotech Major 07)\n\n");
  printf("Hardware by Hector Maciel (Cal Poly Pomona Physics Technician and\n");
  printf("                                           Physics graduate 06)\n\n");
// Get the cpu frequency
  char filename[]="/proc/cpuinfo";
  inf=fopen(filename,"r");
  while(fr<100.0)
   {
     fscanf(inf,"%s", dum);
     fr=atof(dum);
     if(fr>10) frc=fr;
 //    printf("%s    %f\n",dum, fr);
   }

 cpu=frc*1000000.0;
// cpu=1600000000.0;
 printf("The cpu clock frequency is  %f Hz\n\n",cpu);
 printf("menu:\n\n");
    ntot=8;
    ncal=8;
    for (i=0;i<100;i++)
     {
      time1[i]=(i-1)*1.0;
      vel[i]=(i-1)*1.0;
      rvel[i]=vel[i];       
     }
    iopl(3);
    circum=35.0;
    fac=4294967296.0;
    facpic=0.0000004;
    menu();
   running=1;
 do
 {
 
   {
     printf("input an option (then hit return): \n");
     scanf("%c", &cho);
     switch(cho)
     {
      case 'd' : takedat();     //takedata
                 break;
      case 'g' : velplot();    //velplot
                 menu();
                 break;
      case 'v' : datview();    //datview
                 menu;
                 break;
      case 'z' : calbrat();
                 menu;
                 break;
      case 'm' : menu();
                 break;
      case 'q' : running = 2;
      default :  menu();
                 break;
     }


   }
 }while (running==1);
 return(0);
}


void menu()
{
 printf(" total blockings = %d \n", ntot);
 printf(" calibration distance = %f  for %d  blockings\n", circum,ncal);
 printf(" press d to take data\n");
 printf(" press z to change the number of points and calibration distance\n");
 printf(" press g to plot v vs. t\n");
 printf(" press v to view data\n");
 printf(" press m to see the menu\n");
 printf(" press Q to exit the program\n\n\n");
}

void calbrat()
{
 printf("Input the number of blockings \n");
 scanf("%d",&ncal);
 if(ncal<4)ncal=4;
 if(ncal>30)ncal=30;
 ntot=ncal;
 printf("Input the distance between the first and the %d th blocking  \n",ncal);
 scanf("%f",&circum);
 return;
}

void testgate()
 {
 printf(" \n The test is over\n");
 printf("Press t to test again or m to menu\n\n");
 }

void setpoints()
{
  printf("Input the total # of times you want the photogate to be blocked\n");
  printf("for the data analysis. \n");
  scanf("%d",&s1);
  ntot=s1;
  if (s1<ncal) ntot=ncal;
  if (s1>30) ntot=30;
  return;
}

void datview()
 {
 printf("Data Points    Time(sec)         Velocity(cm/sec) \n");
 for (i=1; i<=ntot; i++)
  {
  printf("     %d         %f            %f \n",i, time1[i], rvel[i]);
  }
 }

void velplot()
 {

 #define NIL (0)
 char str1[40];

 Display *dpy = XOpenDisplay(NULL);   //Define Display
 assert(dpy);
 Colormap cmap;                       //To set colors
 XColor color,ignore;
 long fgred, fggreen, fgblue;
 long eventmask = ButtonPressMask|ExposureMask|KeyPressMask|StructureNotifyMask;  //events for the window
 XEvent event;
 int blackColor = BlackPixel(dpy,DefaultScreen(dpy));     //Define colors
 int whiteColor = WhitePixel(dpy,DefaultScreen(dpy));
 cmap = DefaultColormap(dpy,DefaultScreen(dpy));
 XAllocNamedColor(dpy,cmap,"red",&color,&ignore);
 fgred = color.pixel;
 XAllocNamedColor(dpy,cmap,"green",&color,&ignore);
 fggreen = color.pixel;
 XAllocNamedColor(dpy,cmap,"blue",&color,&ignore);
 fgblue = color.pixel;

 Window w = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy),
                        0,0,700,500,0,whiteColor,whiteColor);      //Define window
 XSelectInput(dpy,w,StructureNotifyMask);
 XMapWindow(dpy,w);                                            //Map the window
 GC gc = XCreateGC(dpy,w,0,NIL);                              //Graphics Console
 XSetForeground(dpy,gc,blackColor);

for(;;)
 {
  XEvent e;                                            //Notify that the window is set up
  XNextEvent(dpy,&e);
  if (e.type==MapNotify)
   break;
 }



  scalet = 1;
  scalev = 1;
  maxxx = 0;
  maxxy = 0;
  for (i=1; i<=ntot;i++)
   {
    if (time1[i] > maxxx) maxxx = time1[i];
    if (rvel[i]  > maxxy) maxxy = rvel[i];
   }

  maxx=700;
  maxy=500;

 scalet=0.95*(maxx-200.0)/maxxx;
 delx=maxxx/(maxx-200.0);
 scalev=0.95*(maxy-100.0)/maxxy;
 dely=maxxy/(maxy-100.0);


 sumx = 0.0;
 sumxx = 0.0;
 sumy = 0.0;
 sumxy = 0.0;
 for (i=1; i<=ntot;i++)
     {
        sumx = sumx + time1[i];
        sumxx = sumxx + time1[i]*time1[i];
        sumy = sumy + rvel[i];
        sumxy = sumxy + time1[i]*rvel[i];
     }
 del = ntot*sumxx - sumx*sumx;
 a = (sumxx*sumy-sumx*sumxy)/del;
 b = (ntot*sumxy-sumx*sumy)/del;

 v1c=a+b*time1[1];
 v2c=a+b*time1[ncal];
 denom=v2c*v2c-v1c*v1c;
 if (denom==0) denom=1.0;
   xa=2.0*circum*b/denom;

//begin graphing

  sprintf(str1,"The slope is %10.2f cm/sec^2",xa*b);
  XDrawString(dpy,w,gc,400,450,str1,strlen(str1));
  sprintf(str1,"The intercept is %10.2f cm/sec",xa*a);
  XDrawString(dpy,w,gc,400,480,str1,strlen(str1));

 XDrawLine(dpy,w,gc,200,50,200,400);
 XDrawLine(dpy,w,gc,200,400,700,400);
 char str2[]="Time";
 XDrawString(dpy,w,gc,650,420,str2,strlen(str2));
 char str3[]="Velocity";
 XDrawString(dpy,w,gc,200,40,str3,strlen(str3));
 XFlush(dpy);

 XSetForeground(dpy,gc,fgblue);

 for (i=1; i<=ntot; i++)
  {
    nx=floor(200.0+time1[i]*scalet);
    ny=floor(maxy-100.0-scalev*rvel[i]);
    XDrawArc(dpy,w,gc,nx-5,ny-5,10,10,0,21600);
  }
 XSetForeground(dpy,gc,fggreen);

 for (i=1; i<(maxx-80);i++)
   {
     nx=floor(200+i*delx*scalet);
     ny=floor(maxy-100-scalev*(b*i*delx+a));
     XDrawArc(dpy,w,gc,nx-1,ny-1,2,2,0,21600);
   }

 XSetForeground(dpy,gc,fgred);

 sprintf(str1,"  Time        Velocity");
 XDrawString(dpy,w,gc,5,15,str1,strlen(str1));

 for (i=1;i<=ntot;i++)
 {
  sprintf(str1,"%10.2f  %10.2f",time1[i],rvel[i]);
  XDrawString(dpy,w,gc,5,20+10*i,str1,strlen(str1));
 }

 char str[] = "Press any key to clear graphing window";
 XDrawString(dpy,w,gc,40,(maxy-20),str,strlen(str));



 XSelectInput(dpy,w,eventmask);
 XMapWindow(dpy,w);
 graphing=1;
 do { 
   XWindowEvent(dpy,w,eventmask,&event);
   switch (event.type) {
   case Expose: 
            XSetForeground(dpy,gc,blackColor);
            XDrawLine(dpy,w,gc,200,50,200,400);
            XDrawLine(dpy,w,gc,200,400,700,400);
            char str2[]="Time";
            XDrawString(dpy,w,gc,650,420,str2,strlen(str2));
            char str3[]="Velocity";
            XDrawString(dpy,w,gc,200,40,str3,strlen(str3));
            XFlush(dpy);
            XSetForeground(dpy,gc,fgblue);
            for (i=1; i<=ntot; i++)
             {
                nx=floor(200.0+time1[i]*scalet);
                ny=floor(maxy-100.0-scalev*rvel[i]);
                XDrawArc(dpy,w,gc,nx-5,ny-5,10,10,0,21600);
             }
             XSetForeground(dpy,gc,fggreen);
             for (i=1; i<(maxx-80);i++)
             {
              nx=floor(200+i*delx*scalet);
              ny=floor(maxy-100-scalev*(b*i*delx+a));
              XDrawArc(dpy,w,gc,nx-1,ny-1,2,2,0,21600);
             }
             XSetForeground(dpy,gc,fgred);
              sprintf(str1,"  Time        Velocity");
             XDrawString(dpy,w,gc,5,15,str1,strlen(str1));
             for (i=1;i<=ntot;i++)
             {
              sprintf(str1,"%10.2f  %10.2f",time1[i],rvel[i]);
               XDrawString(dpy,w,gc,5,20+10*i,str1,strlen(str1));
              }
              sprintf(str1,"The slope is %10.2f cm/sec^2",xa*b);
              XDrawString(dpy,w,gc,400,450,str1,strlen(str1));
              sprintf(str1,"The intercept is %10.2f cm/sec",xa*a);
              XDrawString(dpy,w,gc,400,480,str1,strlen(str1));

               char str[] = "Press any key to clear graphing window";
              XDrawString(dpy,w,gc,40,(maxy-20),str,strlen(str));

            break;    
   case KeyPress:  graphing=2;
                   break;       
                                        //Wait till keypressed
  }
 } while(graphing==1);

 XFreeGC(dpy,gc);
 XCloseDisplay(dpy);

 }

void takedat()
{
 asm volatile("rdtsc":"=a"(dum1),"=d"(dum2));
 ncheck=1;

/* open the device to be non-blocking (read will return immediatly) */
        fd = open(MODEMDEVICE, O_RDWR | O_NOCTTY | O_NONBLOCK );
//        if (fd <0) {perror(MODEMDEVICE); exit(-1); }

/* install the signal handler before making the device asynchronous */
        saio.sa_handler = signal_handler_IO;
        sigemptyset(&saio.sa_mask);   //saio.sa_mask = 0;
        saio.sa_flags = 0;
        saio.sa_restorer = NULL;
        sigaction(SIGIO,&saio,NULL);

        /* allow the process to receive SIGIO */
        fcntl(fd, F_SETOWN, getpid());
        /* Make the file descriptor asynchronous (the manual page says only
           O_APPEND and O_NONBLOCK, will work with F_SETFL...) */
        fcntl(fd, F_SETFL, FASYNC);

        tcgetattr(fd,&oldtio); /* save current port settings */
        /* set new port settings for canonical input processing */
        newtio.c_cflag = BAUDRATE | CRTSCTS | CS8 | CLOCAL | CREAD;
        newtio.c_iflag = IGNPAR | ICRNL ;
        newtio.c_oflag = 0;
//        newtio.c_cflag &=~(PARENB | CSTOPB | CSIZE);
//        newtio.c_lflag = ~ICANON;
        newtio.c_cc[VMIN]=8*ntot;
        newtio.c_cc[VTIME]=0;
        newtio.c_lflag &= ~(ICANON | ECHO | ECHOE | ISIG );
        tcflush(fd, TCIFLUSH);
        tcsetattr(fd,TCSANOW,&newtio);
  STOP=FALSE;
  while (STOP==FALSE)
    {
     /* after receiving SIGIO, wait_flag = FALSE, input is available
     and can be read */
     if (wait_flag==FALSE)
          {
            res = read(fd,&buf,255);
            wait_flag = TRUE;      /* wait for new input */
            STOP=TRUE;
          }
      }
     /* restore old port settings */
   asm volatile("rdtsc":"=a"(dum3),"=d"(dum4));
   tcsetattr(fd,TCSANOW,&oldtio);
//  for (i=0;i<8*ntot;i=i+4)
//   printf("%u %u %u %u \n",buf[i],buf[i+1],buf[i+2],buf[i+3]);
  i=0;
  ii=1;
  while(ii<=ntot)
   {
   time1[ii]=buf[i]+256*buf[i+1]+256*256*buf[i+2]+256*256*256*buf[i+3];
   vel[ii]=buf[i+4]+buf[i+5]*256+buf[i+6]*256*256+buf[i+7]*256*256*256;
   vel[ii]=vel[ii]-time1[ii];
   time1[ii]=time1[ii]+vel[ii]/2.0;
   time1[ii]=time1[ii]*facpic;
   vel[ii]=1.0/vel[ii]/facpic;
   ii=ii+1;
   i=i+8;
   }

  tottime=((dum3-dum1)+fac*(dum4-dum2))/cpu;
  toffset=time1[1];
  for(ii=1;ii<=ntot;ii++)
   time1[ii]=time1[ii]-toffset;

  sumx=0.0;
  sumxx=0.0;
  sumy=0.0;
  sumxy=0.0;
  for (i=1; i<=ntot; i++)
  {
  sumx=sumx+time1[i];
  sumxx=sumxx+time1[i]*time1[i];
  sumy=sumy+vel[i];
  sumxy=sumxy+time1[i]*vel[i];
  }
  del=ntot*sumxx-sumx*sumx;
  a=(sumxx*sumy-sumx*sumxy)/del;
  b=(ntot*sumxy-sumx*sumy)/del;

  v1c=a+b*time1[1];
  v2c=a+b*time1[ncal];
  denom=v2c*v2c-v1c*v1c;
  if (denom==0) denom=1;
  xa=2.0*circum*b/denom;

  for (i=1; i<=ntot; i++)
  {
  rvel[i]=xa*vel[i];
  }

  printf("Data Point         Time(sec)         Velocity(cm/sec) \n");
  for (i=1; i<=ntot; i++)
  {
  printf("     %d             %f            %f \n",i, time1[i], xa*vel[i]);
  }
  printf("The acceleration is %f cm/sec2 \n",xa*b);

  {if (runc==2) printf("Bad Data");}
  close(fd);
}

void signal_handler_IO (int status)
 {
   wait_flag = FALSE;
 }