// include the headers of all sequence classes
#include <odinseq/seqall.h>

// The whole EPI sequence is a C++ class
class METHOD_CLASS : public SeqMethod {

private:
  // this section contains the data elements of the sequence,
  // i.e. the sequence specific parameters and the sequence objects

  // the following parameters are specific to EPI sequences,
  // parameters common to all sequences can be accessed via the
  // 'commonPars' global object
  LDRbool   TakeMinEchoTime;
  LDRenum   SliceTimingMode;
  LDRenum   TemplateScan;
  LDRbool   RampSampling;
  LDRenum   RampMode;
  LDRfloat  RampSteepness;
  LDRbool   FatSaturation;
  LDRenum   EPIMode;
  LDRint    Shots;
  LDRint    DummyCycles;
  LDRint    NumOfGradEchoes;
  LDRint    NumOfSamples;
  LDRdouble PulseDur;
  LDRbool   fMRITrigger;
  LDRint    fMRIBlockSize;
  LDRbool   FieldMap;


  // sequence objects which are the elementary objects
  // to build the EPI sequence
  SeqPulsar  exc;
  SeqPulsar  refoc;
  SeqSat     fatsat;
  SeqAcqEPI  epiacq;
  SeqDelay   epiacq_dummy;
  SeqAcqEPI  epiacq_template;
  SeqAcqEPI  epiacq_grappa;
  SeqAcqDeph deph;
  SeqAcqDeph deph_template;
  SeqAcqDeph deph_grappa;
  SeqObjLoop sliceloop;
  SeqObjLoop reploop;
  SeqObjLoop segloop;
  SeqObjLoop grappaloop;
  SeqObjLoop dummyloop;
  SeqDelay   delayuniform;
  SeqDelay   delaybunched;
  SeqObjList scan;
  SeqObjList dummypart;
  SeqObjList templatepart;
  SeqObjList grappapart;
  SeqObjList imagingpart;
  SeqObjList slicepart;
  SeqObjList slicepart_dummy;
  SeqObjList slicepart_template;
  SeqObjList slicepart_grappa;
  SeqObjList preppart;
  SeqDelay   exc2acq;
  SeqDelay   exc2refoc;
  SeqDelay   refoc2acq;
  SeqGradConstPulse spoiler;
  SeqTrigger trigger;
  SeqGradTrapezParallel crusher;
  SeqDelay   crusherdelay;
  SeqFieldMap fmapscan;
  SeqVecIter phaseiter;



public:

// This constructor creates an empty EPI sequence
METHOD_CLASS(const STD_string& label) : SeqMethod(label) {
  set_description("This is a feature-rich multislice EPI method. "
                  "The template scan for phase correction is sampled before the repetitions. "
                  "A field-map scan can be acquired prior to the actual scan for distortion correction.");
}

unsigned int numof_testcases() const {return 3;}

void method_pars_init() {

  // In this function, parameters are initialized and default values are set

  commonPars->set_MatrixSize(readDirection,64);
  commonPars->set_MatrixSize(phaseDirection,64,noedit);
  commonPars->set_NumOfRepetitions(1);
  commonPars->set_RepetitionTime(1000.0);
  commonPars->set_AcqSweepWidth(100.0);

  SliceTimingMode.add_item("EPIuniform");
  SliceTimingMode.add_item("EPIbunchedEarly");
  SliceTimingMode.set_actual("EPIuniform");
  SliceTimingMode.set_description("Scanning of slices uniformly distributed over TR (EPIuniform) or bunched at the start (EPIbunchedEarly)");

  EPIMode.add_item("FIDMode");
  EPIMode.add_item("SEMode");
  EPIMode.set_actual("FIDMode");
  EPIMode.set_description("Gradient-echo EPI (FIDMode) or spin-echo EPI (SEMode)");

  TemplateScan.add_item("NoCorrection");
  TemplateScan.add_item("PhaseCorrection");
  TemplateScan.set_actual("PhaseCorrection");
  TemplateScan.set_description("The type of template scan which is acquired beforehand");

  RampSampling=false;
  RampSampling.set_description("Perform sampling during gradient ramps");

  RampMode.add_item("linear",linear);
  RampMode.add_item("sinusoidal",sinusoidal);
  RampMode.add_item("half_sinusoidal",half_sinusoidal);
  RampMode.set_actual(linear);
  RampMode.set_description("The shape of the ramps of the read gradient");

  RampSteepness=1.0;
  RampSteepness.set_description("Relative steepness (slew rate) of the EPI readout ramps");


  FatSaturation=true;
  FatSaturation.set_description("Saturation of fat resonance prior to excitation");

  TakeMinEchoTime=true;
  TakeMinEchoTime.set_description("Use minimum possible TE");

  Shots=1;
  Shots.set_description("Number of shots, multi-shot (segmented) EPI if Shots>1");

  DummyCycles=3;
  DummyCycles.set_description("Number of dummy shots before actual acquisition");

  fMRITrigger=true;
  fMRITrigger.set_description("External triggering");

  fMRIBlockSize=0;
  fMRIBlockSize.set_description("If non-zero, fMRI analysis will be performed with this size for both stimulation and rest");

  FieldMap=false;
  FieldMap.set_description("Fieldmap pre-scan for distortion correction");

  PulseDur=4.0; // avoid initial high RF amplitude
  PulseDur.set_description("Pulse duration of excitation/refocusing pulse");


  // alternative default settings for sequence test
  if(get_current_testcase()==1) {
    EPIMode.set_actual("SEMode");
    RampSteepness=0.8;
    RampSampling=true;
    FatSaturation=false;
    DummyCycles=0;
    commonPars->set_PartialFourier(0.7);
  }
  if(get_current_testcase()==2) {
    TemplateScan.set_actual("NoCorrection");
    Shots=3;
    FieldMap=true;
    commonPars->set_RFSpoiling(false); 
  }


  // register method parameters for user interface, parameter files, etc.
  append_parameter(EPIMode,"EPIMode");
  append_parameter(TakeMinEchoTime,"TakeMinEchoTime");
  append_parameter(Shots,"Shots");
  append_parameter(DummyCycles,"DummyCycles");
  append_parameter(TemplateScan,"TemplateScan");
  append_parameter(RampSampling,"RampSampling");
  append_parameter(FatSaturation,"FatSaturation");
  append_parameter(FieldMap,"FieldMap");

  fmapscan.init("fmapscan");
  append_parameter(fmapscan.get_parblock(),"FieldMapPars");

  append_parameter(fMRITrigger,"fMRITrigger");
  append_parameter(fMRIBlockSize,"fMRIBlockSize");

  append_parameter(PulseDur,"PulseDur");

  if(systemInfo->get_platform()!=numaris_4) {
    append_parameter(RampMode,"RampMode");
    append_parameter(RampSteepness,"RampSteepness");
    append_parameter(SliceTimingMode,"SliceTimingMode");
    append_parameter(NumOfGradEchoes,"NumOfGradEchoes",noedit);
    append_parameter(NumOfSamples,"NumOfSamples",noedit);
  }


}



void method_seq_init() {

  ///////////////// Pulses: /////////////////////

  float slicethick=geometryInfo->get_sliceThickness();
  float slicegap=geometryInfo->get_sliceDistance()-slicethick;

  // excitation pulse
  exc=SeqPulsarSinc("exc",slicethick,true,PulseDur,commonPars->get_FlipAngle());
  exc.set_rephased(true, 0.8*systemInfo->get_max_grad()); // short rephaser
  exc.set_freqlist( systemInfo->get_gamma() * exc.get_strength() / (2.0*PII) * geometryInfo->get_sliceOffsetVector() );
  exc.set_pulse_type(excitation);

  // Slightly thicker refocusing slice for better SNR
  // Since the same spatial resolution is used for exc and refoc, the gradient strengths will be the same
  float extra_slicethick_refoc=STD_min(0.5*slicethick, 0.3*slicegap);

  // refocusing pulse
  refoc=SeqPulsarSinc("refoc",slicethick+extra_slicethick_refoc,false,PulseDur,180.0);
  refoc.set_freqlist( systemInfo->get_gamma() * refoc.get_strength() / (2.0*PII) * geometryInfo->get_sliceOffsetVector() );
  if(!commonPars->get_RFSpoiling()) refoc.set_phase(90.0);
  refoc.set_pulse_type(refocusing);

  // fat saturation module
  fatsat=SeqSat("fatsat",fat);


  //////////////// EPI-Readout: //////////////////////////////

  // set equivalent resolution in read and phase direction
  float resolution=secureDivision(geometryInfo->get_FOV(readDirection),commonPars->get_MatrixSize(readDirection));
  int pelines=int(secureDivision(geometryInfo->get_FOV(phaseDirection),resolution)+0.5);

  float os_read=1.2; // slight oversampling to allow off-center FOV

  epiacq=SeqAcqEPI("epiacq",commonPars->get_AcqSweepWidth(),
                   commonPars->get_MatrixSize(readDirection),  geometryInfo->get_FOV(readDirection),
                   pelines, geometryInfo->get_FOV(phaseDirection),
                   Shots, commonPars->get_ReductionFactor(), os_read, "", 0, 0, rampType(int(RampMode)),
                   RampSampling, RampSteepness, commonPars->get_PartialFourier());



  // display sampling extents in read/phase direction
  NumOfGradEchoes=epiacq.get_numof_gradechoes();
  commonPars->set_MatrixSize(phaseDirection,NumOfGradEchoes*Shots*commonPars->get_ReductionFactor(),noedit);
  NumOfSamples=epiacq.get_npts_read();

  // Template scan fo EPI, begin with copy of actual EPI
  epiacq_template=epiacq;
  epiacq_template.set_label("epiacq_template");

  // 1D phase correction
  if(TemplateScan=="PhaseCorrection") {
    epiacq_template.set_template_type(phasecorr_template);
  }


  // Full multi-shot EPI readout as GRAPPA training data
  epiacq_grappa=epiacq;
  epiacq_grappa.set_label("epiacq_grappa");
  epiacq_grappa.set_template_type(grappa_template);


  // Delay instead of actual EPI readout for dummy scans
  epiacq_dummy=SeqDelay("epiacq_dummy",epiacq.get_duration());


  // EPI pre-dephase gradient
  dephaseMode dephmode=FID;
  if (EPIMode=="SEMode") dephmode=spinEcho;
  deph=SeqAcqDeph("deph",epiacq,dephmode);
  deph_template=SeqAcqDeph("deph_template",epiacq_template,dephmode);
  deph_grappa=SeqAcqDeph("deph_grappa",epiacq_grappa,dephmode);



  /////////////////// RF Spoiling ///////////////////////////////////////////////////////

  if(commonPars->get_RFSpoiling()) {

     // recommended by Goerke et al., NMR Biomed. 18, 534-542 (2005)
    int plistsize=16;
    double plistincr=45.0;

    exc.set_phasespoiling(plistsize, plistincr);
    refoc.set_phasespoiling(plistsize, plistincr, 90.0);

    epiacq.set_phasespoiling(plistsize, plistincr);
    epiacq_template.set_phasespoiling(plistsize, plistincr);
    epiacq_grappa.set_phasespoiling(plistsize, plistincr);

    phaseiter=SeqVecIter("phaseiter");
    phaseiter.add_vector(exc.get_phaselist_vector());
    phaseiter.add_vector(refoc.get_phaselist_vector());
    phaseiter.add_vector(epiacq.get_phaselist_vector());
    phaseiter.add_vector(epiacq_template.get_phaselist_vector());
    phaseiter.add_vector(epiacq_grappa.get_phaselist_vector());
  }

  //////////////// Loops: //////////////////////////////

  // loop to iterate over slices
  sliceloop=SeqObjLoop("sliceloop");

  // loop to iterate over repetitions
  reploop=SeqObjLoop("reploop");

  // loop to iterate over segments (multi shot)
  segloop=SeqObjLoop("segloop");

  // loop to iterate over GRAPPA interleaves to obtain training data
  grappaloop=SeqObjLoop("grappaloop");

  // loop to iterate over dummy scans
  dummyloop=SeqObjLoop("dummyloop");

  //////////////// Timing Delays: //////////////////////////////

  // Delays uniform/bunched
  delayuniform=SeqDelay("delayuniform");
  delaybunched=SeqDelay("delaybunched");


  //////////////// Spoiler Gradient: //////////////////////////////

  double spoiler_strength=0.5*systemInfo->get_max_grad();
  double spoiler_integral=4.0*fabs(deph.get_gradintegral().sum());
  double spoiler_dur=secureDivision(spoiler_integral,spoiler_strength);

  spoiler=SeqGradConstPulse("spoiler",readDirection,spoiler_strength,spoiler_dur);


  //////////////// Crusher Gradient: //////////////////////////////

  float crusher_integral=2.0*spoiler_integral;
  crusher=SeqGradTrapezParallel("crusher",crusher_integral,crusher_integral,crusher_integral, spoiler_strength);

  crusherdelay=SeqDelay("crusherdelay",0.1); // Small delay to avoid gradient-induced stimulation


  //////////////// fMRI trigger: //////////////////////////////

  trigger=SeqTrigger("fmri_trigger",1.0);


  //////////////// Field-map template: //////////////////////////////

  if(FieldMap) {
    if(FatSaturation) fmapscan.build_seq(commonPars->get_AcqSweepWidth(),os_read,fatsat); // pass fat saturation on to field-map scan
    else              fmapscan.build_seq(commonPars->get_AcqSweepWidth(),os_read);
  }


  //////////////// Build the sequence: //////////////////////////////


  // add fat saturation to template and repetitions
  if(FatSaturation) preppart += fatsat;

  if(fMRITrigger) slicepart += trigger;

  if (EPIMode=="FIDMode") {

    dummypart=    preppart + exc + deph          + exc2acq + epiacq_dummy     + crusherdelay + crusher;

    templatepart= preppart + exc + deph_template + exc2acq + epiacq_template  + crusherdelay + crusher;

    grappapart=   preppart + exc + deph_grappa   + exc2acq + epiacq_grappa    + crusherdelay + crusher;

    imagingpart=  preppart + exc + deph          + exc2acq + epiacq           + crusherdelay + crusher;


    slicepart_dummy =    sliceloop( dummypart    + delayuniform )[exc];

    slicepart_template = sliceloop( templatepart + delayuniform )[exc];

    slicepart_grappa =   sliceloop( grappapart   + delayuniform )[exc];

    slicepart +=         sliceloop( imagingpart  + delayuniform )[exc];
  }


  if (EPIMode=="SEMode") {

    dummypart=    preppart + exc + exc2refoc + deph          + spoiler + refoc + spoiler + refoc2acq + epiacq_dummy    + crusherdelay + crusher;

    templatepart= preppart + exc + exc2refoc + deph_template + spoiler + refoc + spoiler + refoc2acq + epiacq_template + crusherdelay + crusher;

    grappapart=   preppart + exc + exc2refoc + deph_grappa   + spoiler + refoc + spoiler + refoc2acq + epiacq_grappa   + crusherdelay + crusher;

    imagingpart=  preppart + exc + exc2refoc + deph          + spoiler + refoc + spoiler + refoc2acq + epiacq          + crusherdelay + crusher;


    slicepart_dummy =    sliceloop( dummypart + delayuniform    )[exc][refoc];

    slicepart_template = sliceloop( templatepart + delayuniform )[exc][refoc];

    slicepart_grappa  =  sliceloop( grappapart + delayuniform   )[exc][refoc];

    slicepart +=         sliceloop( imagingpart + delayuniform  )[exc][refoc];
  }

  if(commonPars->get_RFSpoiling()) {
    slicepart          += phaseiter;
    slicepart_dummy    += phaseiter;
    slicepart_template += phaseiter;
    slicepart_grappa   += phaseiter;
  }

  if(FieldMap) scan += fmapscan + delayuniform + delaybunched;


  if(DummyCycles>0) {
    scan+= dummyloop(
             slicepart_dummy + delaybunched
           )[DummyCycles];
  }



  if(TemplateScan=="PhaseCorrection") {
    scan += slicepart_template + delaybunched; // phasecorr requires only one shot from template
  }


  if(commonPars->get_ReductionFactor()>1) {
    // Fully sampled k-space
    scan+= grappaloop(
             segloop(
               slicepart_grappa + delaybunched
             )[deph_grappa.get_epi_segment_vector()]
           )[deph_grappa.get_epi_reduction_vector()];
  }


  scan+= reploop(
           segloop(
             slicepart + delaybunched
           )[deph.get_epi_segment_vector()]
         )[commonPars->get_NumOfRepetitions()];


  set_sequence( scan );
}




void method_rels() {


   ////////////////// TE Timings: ////////////////////////////////


   double min_echo_time=0.0;


   if (EPIMode=="FIDMode") {

     min_echo_time=(exc.get_duration()-exc.get_magnetic_center())+deph.get_duration()+epiacq.get_acquisition_center();

     if(commonPars->get_EchoTime()<min_echo_time) commonPars->set_EchoTime(min_echo_time);

     if(TakeMinEchoTime) commonPars->set_EchoTime(min_echo_time);

     exc2acq=commonPars->get_EchoTime()-min_echo_time;
   }



   if (EPIMode=="SEMode") {

     //////////////// Duration from the middle of the excitation pulse ///////////////////
     //////////////// to the middle of the refocusing pulse: /////////////////////////////

     double TE1=( exc.get_duration() - exc.get_magnetic_center() )
                + deph.get_duration()
                + spoiler.get_duration()
                + refoc.get_magnetic_center();

     if(min_echo_time<(2.0*TE1)) min_echo_time=2.0*TE1;

     //////////////// Duration from the middle of the refocusing pulse ///////////////////
     //////////////// to the middle of the acquisition window: ///////////////////////////

     double TE2=(refoc.get_duration() - refoc.get_magnetic_center())
                + spoiler.get_duration()
                + epiacq.get_acquisition_center();

     if(min_echo_time<(2.0*TE2)) min_echo_time=2.0*TE2;



     if (commonPars->get_EchoTime()<min_echo_time) commonPars->set_EchoTime(min_echo_time);

     if(TakeMinEchoTime) commonPars->set_EchoTime(min_echo_time);

     exc2refoc=commonPars->get_EchoTime()/2.0-TE1;
     refoc2acq=commonPars->get_EchoTime()/2.0-TE2;

   }



   ////////////////// TR Timings: ////////////////////////////////

   // reset before recalculating timings
   delayuniform=0.0;
   delaybunched=0.0;

   float slicedur=slicepart.get_duration();
   if(commonPars->get_RepetitionTime()<slicedur) commonPars->set_RepetitionTime(slicedur);

   if (SliceTimingMode=="EPIuniform") {
     delayuniform=(commonPars->get_RepetitionTime()-slicedur)/double(geometryInfo->get_nSlices());
     delaybunched=0.0;
   }

   if (SliceTimingMode=="EPIbunchedEarly") {
     delaybunched=commonPars->get_RepetitionTime()-slicedur;
     delayuniform=0.0;
   }


}


void method_pars_set() {

  // extra information for the automatic reconstruction
  epiacq.set_reco_vector(slice,exc);
  epiacq_template.set_reco_vector(slice,exc);
  epiacq_grappa.set_reco_vector(slice,exc);

  if(fMRIBlockSize>0) recoInfo->set_PostProc3D("repcoll | fmri("+itos(fMRIBlockSize)+","+itos(fMRIBlockSize)+")");

}

};


/////////////////////////////////////////////////////


// entry point for the sequence module
ODINMETHOD_ENTRY_POINT