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


// The whole EPI sequence (including reco) is a C++ class
class METHOD_CLASS : public SeqMethod {

private:

  LDRint   NumOfEchoes;
  LDRfloat PulseDur;
  LDRint   NumOfSatPulses;

  LDRdoubleArr TEs;

  // sequence objects which are the elementary objects
  // to build the EPI sequence
  SeqPulsar  exc;
  SeqPulsarReph exc_reph;
  SeqPulsar  refoc;
  SeqSat     fatsat;

  SeqAcqRead  ge_acq;
  SeqAcqRead  se_acq;

  SeqAcqDeph deph;
  SeqObjList ge_preacq;
  SeqObjList se_preacq;

  SeqAcqDeph reph;

  SeqDelay   sedelay;

  SeqGradTrapezParallel rewind;

  SeqGradPhaseEnc pe,pe_rewind;

  SeqVecIter phaseiter;

  SeqObjLoop sliceloop;
  SeqObjLoop peloop;
  SeqObjLoop reploop;
  SeqObjLoop echoloop;
  SeqObjLoop dummyloop;

  SeqDelay   trdelay;

  SeqObjVector gerewindvec;
  SeqObjList gepart;
  SeqDelay   gepart_dummy;

  SeqObjVector serewindvec;
  SeqObjList separt;
  SeqDelay   separt_dummy;

  SeqObjList scan;
  SeqObjList imagingpart;
  SeqObjList dummypart;
  SeqObjList slicepart;
  SeqObjList slicepart_dummy;
  SeqObjList preppart;


  SeqDelay   exc2acq;
  SeqDelay   exc2refoc;
  SeqDelay   refoc2acq;
  SeqGradConstPulse spoiler;

  SeqGradTrapezParallel crusher;
  SeqDelay   crusherdelay;


public:

// This constructor creates an empty EPI sequence
METHOD_CLASS(const STD_string& label) : SeqMethod(label) {
  set_description("GESSE sequence.");
}


void method_pars_init() {

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

  commonPars->set_MatrixSize(readDirection,128);
  commonPars->set_MatrixSize(phaseDirection,128,noedit);
  commonPars->set_NumOfRepetitions(1);
  commonPars->set_RepetitionTime(1000.0);
  commonPars->set_AcqSweepWidth(100.0);
//  commonPars->set_ReductionFactor(3);

  NumOfEchoes=32;
  NumOfEchoes.set_description("Number of echoes per period");
  append_parameter(NumOfEchoes,"NumOfEchoes");

  PulseDur=4.0; // start with large duration so that sequence can be loaded
  PulseDur.set_unit(ODIN_TIME_UNIT);
  PulseDur.set_description("Pulse duration of excitation and refocusing pulse");
  append_parameter(PulseDur,"PulseDur");

  NumOfSatPulses=2;
  NumOfSatPulses.set_description("Number of consecutive saturation pulses");
  append_parameter(NumOfSatPulses,"NumOfSatPulses");

  TEs.resize(2*NumOfEchoes);
  TEs.set_description("Echo times");
  append_parameter(TEs,"TEs");

}



void method_seq_init() {
  Log<Seq> odinlog(this,"method_seq_init");

  if(NumOfEchoes<2) NumOfEchoes=2;

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

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

  // excitation pulse
  exc=SeqPulsarSinc("exc",slicethick,false,PulseDur,commonPars->get_FlipAngle());
  dvector exclist=systemInfo->get_gamma() * exc.get_strength() / (2.0*PII) * geometryInfo->get_sliceOffsetVector();
  ODINLOG(odinlog,normalDebug) << "exclist=" << exclist << STD_endl;
  exc.set_freqlist( exclist );
  exc.set_pulse_type(excitation);

  exc_reph=SeqPulsarReph("exc_reph",exc);

  // 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);
  dvector refoclist=systemInfo->get_gamma() * refoc.get_strength() / (2.0*PII) * geometryInfo->get_sliceOffsetVector();
  ODINLOG(odinlog,normalDebug) << "refoclist=" << refoclist << STD_endl;
  refoc.set_freqlist( refoclist );
  if(!commonPars->get_RFSpoiling()) refoc.set_phase(90.0);
  refoc.set_pulse_type(refocusing);

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



  //////////////// 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);
  commonPars->set_MatrixSize(phaseDirection, pelines, noedit);

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


  ge_acq=SeqAcqRead("ge_acq",commonPars->get_AcqSweepWidth(),commonPars->get_MatrixSize(readDirection),
                               geometryInfo->get_FOV(readDirection),readDirection,os_read);


  se_acq=ge_acq;

  // pre-dephase gradient
  deph=SeqAcqDeph("deph",ge_acq,FID);

  // post-rephase gradients
  reph=SeqAcqDeph("reph",ge_acq,rephase);

  // EPI rewinder after each echo
  fvector gradint(3);

  float rewind_strength=0.4*systemInfo->get_max_grad(); // OK for stimulation monitor 

  gradint=reph.get_gradintegral()+deph.get_gradintegral(); // collapse rephase and dephase into one gradient pulse
  rewind=SeqGradTrapezParallel("rewind",gradint[0],gradint[1],gradint[2],rewind_strength);

  ODINLOG(odinlog,significantDebug) << "rewind.get_duration()=" << rewind.get_duration() << STD_endl;


  //////////////// Phase Encoding: //////////////////////////

  pe=SeqGradPhaseEnc("pe",commonPars->get_MatrixSize(phaseDirection),geometryInfo->get_FOV(phaseDirection),
                  phaseDirection,0.25*systemInfo->get_max_grad(),
                  linearEncoding,noReorder,1,
                  commonPars->get_ReductionFactor(), DEFAULT_ACL_BANDS, commonPars->get_PartialFourier());

  pe_rewind=pe;
  pe_rewind.set_label("pe_rewind");
  pe_rewind.invert_strength();



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

  if(commonPars->get_RFSpoiling()) {

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

    // Handbook of MRI
    int plistsize=80;
    double plistincr=117.0;

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

    ge_acq.set_phasespoiling(plistsize, plistincr);
    se_acq.set_phasespoiling(plistsize, plistincr);


    phaseiter=SeqVecIter("phaseiter");
    phaseiter.add_vector(exc.get_phaselist_vector());
    phaseiter.add_vector(refoc.get_phaselist_vector());
    phaseiter.add_vector(ge_acq.get_phaselist_vector());
    phaseiter.add_vector(se_acq.get_phaselist_vector());

  }


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

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

  // loop to iterate over phase encoding
  peloop=SeqObjLoop("peloop");

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

  // loop to iterate over EPI modules
  echoloop=SeqObjLoop("echoloop");

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


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

  // TR delay
  trdelay=SeqDelay("trdelay");

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

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

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

//  zoomspoiler=SeqGradConstPulse("zoomspoiler",readDirection,spoiler_strength,0.5*spoiler_dur);

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

  float crusher_strength=0.3*systemInfo->get_max_grad(); // Moderate strength to avoid problems with stimulation
  float crusher_integral=2.0*spoiler_integral;
  crusher=SeqGradTrapezParallel("crusher",crusher_integral,crusher_integral,crusher_integral, crusher_strength);

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


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


  // GE part
  ivector ge_iv(NumOfEchoes);
  int i;
  for(i=0; i<NumOfEchoes; i++) {
    if(i<(NumOfEchoes-1)) {
      gerewindvec          += rewind;
    } else {
      gerewindvec          += reph / pe_rewind  / spoiler;
    }
    ge_iv[i]=i; // Assign index 0-(NumOfEchoes-1) to GE part
  }
  gerewindvec.set_indexvec(ge_iv);


  // SE part
  int nechoes_se=2*NumOfEchoes-1; // sample whole spin echo
  ivector se_iv(nechoes_se);
  for(i=0; i<(nechoes_se); i++) {
    if(i<(nechoes_se-1)) {
      serewindvec          += rewind;
    } else {
      serewindvec          += reph / pe_rewind;
    }
    se_iv[i]=NumOfEchoes+i; // Assign index NumOfEchoes-(3*NumOfEchoes-1) to SE part
  }
  serewindvec.set_indexvec(se_iv);

  ge_preacq.set_label("ge_preacq");
  se_preacq.set_label("se_preacq");

  ge_preacq+=deph/pe/exc_reph;
  se_preacq+=deph/pe;

  // Balanced multi-echo trains
  gepart = ge_preacq + echoloop( ge_acq + gerewindvec )[gerewindvec];

  gepart_dummy = SeqDelay("gepart_dummy", gepart.get_duration());

  separt = se_preacq + echoloop( se_acq + serewindvec )[serewindvec];

  separt_dummy = SeqDelay("separt_dummy", separt.get_duration());


  sedelay=SeqDelay("sedelay",0.0);

  imagingpart= preppart + exc + gepart       + refoc + spoiler + sedelay + separt       + crusherdelay + crusher + crusherdelay;

  dummypart=   preppart + exc + gepart_dummy + refoc + spoiler + sedelay + separt_dummy + crusherdelay + crusher + crusherdelay;



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

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


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


  scan += dummyloop( slicepart_dummy )[3];


  // actual scan loop
  scan+= reploop(
           peloop(
             slicepart
           )[pe][pe_rewind]
         )[commonPars->get_NumOfRepetitions()];


  set_sequence( scan );
}




void method_rels() {
  Log<Seq> odinlog(this,"method_rels");


  double TEexc= exc.get_duration() - exc.get_magnetic_center();
  ODINLOG(odinlog,significantDebug) << "TEexc=" << TEexc << STD_endl;


  // Fixed TE according to GE part
  double TEhalf = TEexc
             + gepart.get_duration()
             + refoc.get_magnetic_center();

  commonPars->set_EchoTime(2.0*TEhalf);

  int nechoes_se=2*NumOfEchoes-1; // sample whole spin echo

  TEs.resize(NumOfEchoes+nechoes_se);

  double TEpreacq_ge=ge_preacq.get_duration();
  double TEpreacq_se=se_preacq.get_duration();
  double TEread=ge_acq.get_duration()+rewind.get_duration();
  double TEacq=ge_acq.get_acquisition_center();
  ODINLOG(odinlog,significantDebug) << "TEpreacq_ge/TEpreacq_se/TEread/TEacq=" << TEpreacq_ge << "/" << TEpreacq_se << "/" << TEread << "/" << TEacq << STD_endl;

  // GE echo times
  int i;
  for(i=0; i<NumOfEchoes; i++) {
    TEs[i]=TEexc+TEpreacq_ge+float(i)*TEread+TEacq;
    ODINLOG(odinlog,significantDebug) << "TEs[" << i << "]=" << TEs[i] << STD_endl;
  }

  double TErefoc=refoc.get_duration() - refoc.get_magnetic_center() + spoiler.get_duration();

  // SE echo times  
  for(i=0; i<nechoes_se; i++) {
    TEs[NumOfEchoes+i]=TEhalf+TErefoc+TEpreacq_se+float(i)*TEread+TEacq;
  }


  // Check placement of SE-rephasing echo times
  sedelay=0.0; // reset
  for(i=0; i<NumOfEchoes; i++) {
    if(TEs[NumOfEchoes+i]>commonPars->get_EchoTime()) {
      ODINLOG(odinlog,errorLog) << "TE[" << (NumOfEchoes+i) << "] too late by " << fabs(TEs[NumOfEchoes+i]-commonPars->get_EchoTime()) << STD_endl;
    }
  }

  // middle readout samples SE
  double tediff=commonPars->get_EchoTime()-TEs[2*NumOfEchoes-1];
  sedelay=tediff;
  ODINLOG(odinlog,normalDebug) << "tediff=" << tediff << STD_endl;
  for(i=0; i<nechoes_se; i++) TEs[NumOfEchoes+i]+=tediff;


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

  float slicedur=slicepart.get_duration();
  if(commonPars->get_RepetitionTime()<slicedur) commonPars->set_RepetitionTime(slicedur);
  trdelay=(commonPars->get_RepetitionTime()-slicedur)/double(geometryInfo->get_nSlices());

}


void method_pars_set() {

  // extra information for the automatic reconstruction
  ge_acq.set_reco_vector(slice,exc);
  se_acq.set_reco_vector(slice,exc);
  ge_acq.set_reco_vector(line,pe);
  se_acq.set_reco_vector(line,pe);


  // Index for TEs
  ge_acq.set_reco_vector(userdef,gerewindvec);
  se_acq.set_reco_vector(userdef,serewindvec);
  recoInfo->set_DimValues(userdef,TEs);


  recoInfo->set_PostProc3D("usercoll | messer");

  recoInfo->set_CmdLineOpts("-ff Hamming -fp 0.8");

}

};

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


// entry point for the sequence module
ODINMETHOD_ENTRY_POINT