#include "find_MAT_akc.h"
#include "get_tac_akc.h"
#include "struct_ROI.h"
#include "struct_TAC_in.h"
#include "constants.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/* C source code name: get_tac_akc.c */
/* edit, e.g., with emacs. compile with Makefile  */

/* By: Alexander K. Converse, University of Wisconsin-Madison */
/* Created: 1 Dec 2012 */
/* Modified: 26 March 2014 */

/* Purpose: find Mean Arrival Time for one ROI  */
/* Note: the model is hard coded to use a specified number (n60Frames from input file) */
/* of 60 second frames at the beginning of the scan. The rest of the frames are assumed */
/* to be 300 second frames. */
/* usage: find_MAT_akc(struct ROI, struct TAC, double sv, int MAT) */

void find_MAT_akc(int index, int flag, struct ROI *iROI, struct TAC_t *q, double sv, int MAT)
{
  int i = 0;
  if (flag == 0) {
    printf("Area: %f\n", q->area[index] * 100);
    printf("%s\n", q->ROIs[index]);
    printf("%f\n", q->admAct);
  }

  for (i = 0; i < q->nFrames; ++i) {
    // Set correct frame durations
    if (i < q->n60Frames) {
      iROI->duration[i] = 60.0;
    } else { 
      iROI->duration[i] = 300.0;
    }

    //((Actual * conversionFactor * area * conversionfactor * sliceThickness) / sliceThickness) * administeredActivity
    iROI->observed[i] = ((((iROI->actualObserved[i] / 1000) * q->area[index] * 100 * q->sliceThickness)/ q->sliceThickness)/ q->admAct);

    // Set observed to an adjusted value
    // Observed = Input Observed / Administered Amount
    if (flag == 0) {
      printf("actualObserved: %f\n", iROI->actualObserved[i]);
    }

    // Set 0th bound value to 0; set other values as determined by equation
    // Bound[i] = Bound[i-1]+ Duration[i-1]*sv*(Observed[i-1]-Bound[i-1]) 
    if (i == 0) {
      iROI->bound[i] = 0;
    } else {
      iROI->bound[i] = iROI->bound[i-1] + (iROI->duration[i-1] * sv * (iROI->observed[i-1] - iROI->bound[i-1]));
    }

    // free = observed - bound
    iROI->free[i] = iROI->observed[i] - iROI->bound[i];
    if (iROI->free[i] < 0) {
      iROI->free[i] = 0;
    }

    // M = bound + free
    iROI->M[i] = iROI->bound[i] + iROI->free[i];

    // Set 0th sumG to different value than other sumG values
    // sumG[i] = sumG[i-1] + Duration[i]*Free[i]
    if (i == 0) {
      iROI->sumG[i] = (60) * iROI->free[i];
    } else {
      iROI->sumG[i] = iROI->sumG[i-1] + (iROI->duration[i] * iROI->free[i]);
    }

    // Set 0th sumG*t to different value than other sumG*t values
    // sumGt = sumGt[i-1] + Duration[i]*Free[i]*Time[i]
    if (i == 0) {
      iROI->sumGt[i] = (60) * iROI->free[i] * iROI->time[i];
    } else {
      iROI->sumGt[i] = iROI->sumGt[i-1] + (iROI->duration[i] * iROI->free[i] * iROI->time[i]);
    }
    
  }

  // Calulate MAT
  // MAT = sumGt[MAT] / sumG[MAT]
  iROI->MAT = (iROI->sumGt[MAT] / iROI->sumG[MAT]);
}