// Function to read rawdata into Matlab
//
// Use 'make' on console to generate RawDataCTX.o.
// Use 'mex matlabread.cc RawDataCTX.o' in Matlab to compile.
// 
// When called with an empy filename, the previously open data file will be used.
//
// Oliver Grimm, September 2009


#include <unistd.h>
#include <string.h>

#include "/ihp/local/Matlab08/extern/include/mex.h"
#include "../../drsdaq/RawDataCTX.h"


// static declaration keep variable alive over invocations of mex file
static RawDataCTX *RD = NULL;

// Exit function called when Matlab is terminated (RD==NULL is OK for delete)
void ExitFcn() {
  delete RD;
}

// Interface function to Matlab
void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

  mexAtExit(*ExitFcn); // Exit function to call when matlab terminates
  
  if (RD == NULL) RD = new RawDataCTX(true);	// Instantiate without console output

  // Check inputs and output arguments
  if (nrhs<1 || nrhs>2 || (nrhs==2 && !mxIsDouble(prhs[1])) || nlhs<1 || nlhs>5) {
    mexErrMsgTxt("Usage: [Data TrigCells RHeader BStructs EHeader] = matlabread('Filename', EventNo)");
  }
  
  char *Filename = mxArrayToString(prhs[0]);
  if (Filename == NULL) mexErrMsgTxt("Could not create filename string.");
  
  if (strlen(Filename) != 0) {
    switch (RD->OpenDataFile(Filename), NULL) {
      case CTX_FOPEN:   mexErrMsgTxt("Could not open file.");
      case CTX_RHEADER: mexErrMsgTxt("Could not read run header.");
      case CTX_BSTRUCT: mexErrMsgTxt("Could not read board structures.");
    }
    if (RD->RHeader->MagicNum == MAGICNUM_OPEN) {
      mexWarnMsgTxt("Magic number in run header indicates that the file has not "
    		    "been closed properly.");
    }
    if (RD->RHeader->MagicNum == MAGICNUM_ERROR) {
      mexWarnMsgTxt("Magic number in run header indicates that an error occurred "
    		    "while writing the file.");
    }
  }
  mxFree(Filename);
  
  // ...some abbrevation for convenience
  unsigned int Boards = RD->RHeader->NBoards;
  unsigned int Channels = RD->RHeader->NChannels;
  unsigned int Chips = RD->RHeader->NChips;
  unsigned int Samples = RD->RHeader->Samples;
  
  // Read event
  if(RD->ReadEvent(mxIsDouble(prhs[1]) ? (int) mxGetScalar(prhs[1]):0, NULL) != CTX_OK) {
    mexErrMsgTxt("Could not read event.");
  }

  // Dummy loop to allow different return values (break will be issued depending on nlhs)
  while(true) {

  // ========= Data array ==========
  
  plhs[0] = mxCreateDoubleMatrix(Boards*Chips*Channels*Samples, 1, mxREAL);
  unsigned int Count = 0;
  for(unsigned int i=0; i<Boards; i++) {
    for(unsigned int j=0; j<Chips*Channels*Samples; j++) {
      *(mxGetPr(plhs[0])+Count) = (double) *((short *) (RD->Data+Boards*Chips*sizeof(int)) + Count++) * RD->BStruct[i].ScaleFactor; 
    }
  }
  // Resize into boards, channels and samples array
  mwSize Dimensions[] = {Samples, Chips*Channels, Boards};
  mxSetDimensions(plhs[0], Dimensions, sizeof(Dimensions)/sizeof(mwSize));

  if (nlhs==1) break;

  // ========= Trigger cell array ==========
    
  plhs[1] = mxCreateDoubleMatrix(Chips, Boards, mxREAL);
  for(unsigned int i=0; i<Boards*Chips; i++) {
    *(mxGetPr(plhs[1])+i) = *((int *) (RD->Data+i*sizeof(int))); 
  }
  if (nlhs==2) break;
 
  // ========= Run header structure ==========
    
  const char *RH_Fields[] = {"MagicNum", "DataFormat", "RunHeaderSize", "EventHeaderSize", "BoardStructureSize", "SoftwareRevision", "Identification", "Type","RunNumber", "FileNumber", "Description", "NBoards", "NChips", "NChannels", "Samples", "Offset", "Events", "NBytes", "StartSecond", "StartMicrosecond", "EndSecond", "EndMicrosecond"};
  const mwSize RH_Dim[] = {1, 1};
  plhs[2] = mxCreateStructArray(sizeof(RH_Dim)/sizeof(mwSize), RH_Dim, sizeof(RH_Fields)/sizeof(char *), RH_Fields);
  mxArray *Data[sizeof(RH_Fields)/sizeof(char *)]; 

  Data[0] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[0]), &RD->RHeader->MagicNum, 4);
  Data[1] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[1]), &RD->RHeader->DataFormat, 4);
  Data[2] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[2]), &RD->RHeader->RunHeaderSize, 4);
  Data[3] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[3]), &RD->RHeader->EventHeaderSize, 4);
  Data[4] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[4]), &RD->RHeader->BoardStructureSize, 4);
  Data[5] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[5]), &RD->RHeader->SoftwareRevision, 4);
  Data[6] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[6]), &RD->RHeader->Identification, 4);
  Data[7] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[7]), &RD->RHeader->Type, 4);
  Data[8] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[8]), &RD->RHeader->RunNumber, 4);
  Data[9] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[9]), &RD->RHeader->FileNumber, 4);
  Data[10] = mxCreateString(RD->RHeader->Description);
  Data[11] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[11]), &RD->RHeader->NBoards, 4);
  Data[12] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[12]), &RD->RHeader->NChips, 4);
  Data[13] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[13]), &RD->RHeader->NChannels, 4);
  Data[14] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[14]), &RD->RHeader->Samples, 4);
  Data[15] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[15]), &RD->RHeader->Offset, 4);
  Data[16] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[16]), &RD->RHeader->Events, 4);
  Data[17] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[17]), &RD->RHeader->NBytes, 4);
  Data[18] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[18]), &RD->RHeader->StartSecond, 4);
  Data[19] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[19]), &RD->RHeader->StartMicrosecond, 4);
  Data[20] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[20]), &RD->RHeader->EndSecond, 4);
  Data[21] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  memcpy(mxGetPr(Data[21]), &RD->RHeader->EndMicrosecond, 4);

  for (int i=0; i<sizeof(RH_Fields)/sizeof(char *); i++) {
    mxSetField(plhs[2], 0, RH_Fields[i], Data[i]);
  }

  if (nlhs==3) break;

  // ========= Board structures ==========
    
  const char *BS_Fields[] = {"SerialNo","NomFreq","BoardTemp","ScaleFactor"};
  const mwSize BS_Dim[] = {1, Boards};
  plhs[3] = mxCreateStructArray(sizeof(BS_Dim)/sizeof(mwSize), BS_Dim, sizeof(BS_Fields)/sizeof(char *), BS_Fields);
  mxArray* (*BSData)[sizeof(BS_Fields)/sizeof(char *)] = new mxArray* [Boards][sizeof(BS_Fields)/sizeof(char *)]; 
 
  for (int j=0; j<Boards; j++) {
    BSData[j][0] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
    BSData[j][1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);
    BSData[j][2] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);
    BSData[j][3] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);

    memcpy(mxGetPr(BSData[j][0]), &RD->BStruct[j].SerialNo, 4);
    memcpy(mxGetPr(BSData[j][1]), &RD->BStruct[j].NomFreq, 4);
    memcpy(mxGetPr(BSData[j][2]), &RD->BStruct[j].BoardTemp, 4);
    memcpy(mxGetPr(BSData[j][3]), &RD->BStruct[j].ScaleFactor, 4);

    for (int i=0; i<sizeof(BS_Fields)/sizeof(char *); i++) {
      mxSetField(plhs[3], j, BS_Fields[i], BSData[j][i]);
    }
  }
  if (nlhs==4) break;

  // ========= Event header structure ==========
    
  const char *EH_Fields[] = {"EventNumber", "Second", "Microsecond", "TriggerType", "EventSize"};
  const mwSize EH_Dim[] = {1, 1};
  plhs[4] = mxCreateStructArray(sizeof(EH_Dim)/sizeof(mwSize), EH_Dim, sizeof(EH_Fields)/sizeof(char *), EH_Fields);
 mxArray *EHData[sizeof(EH_Fields)/sizeof(char *)]; 
 
  for (int i=0; i<sizeof(EH_Fields)/sizeof(char *); i++) {
    EHData[i] = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
  }
  memcpy(mxGetPr(EHData[0]), &RD->EHeader->EventNumber, 4);
  memcpy(mxGetPr(EHData[1]), &RD->EHeader->Second, 4);
  memcpy(mxGetPr(EHData[2]), &RD->EHeader->Microsecond, 4);
  memcpy(mxGetPr(EHData[3]), &RD->EHeader->TriggerType, 4);
  memcpy(mxGetPr(EHData[4]), &RD->EHeader->EventSize, 4);
  for (int i=0; i<sizeof(EH_Fields)/sizeof(char *); i++) {
    mxSetField(plhs[4], 0, EH_Fields[i], EHData[i]);
  }
  if (nlhs>=5) break;

  } // Dummy while loop  
}