#include <odinseq/seqall.h>

class METHOD_CLASS : public SeqMethod {

private:
  LDRint    Blades;
  LDRbool   ShortAxis;
  LDRfloat  BladeOversampling;
  LDRbool   TakeMinEchoTime;
  LDRenum   TemplateScan;
  LDRbool   RampSampling;
  LDRenum   RampMode;
  LDRfloat  RampSteepness;
  LDRbool   FatSaturation;
  LDRfloat  T1Ernst;
  LDRenum   EPIMode;
  LDRint    DummyCycles;
  LDRint    NumOfGradEchoes;
  LDRint    NumOfSamples;
  LDRdouble PulseDur;
  LDRbool   fMRITrigger;
  LDRbool   UpDownBlade;
  LDRbool   FieldMap;
  LDRbool   InvertPartialFourier;


  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 bladeloop;
  SeqObjLoop grappaloop;
  SeqObjLoop dummyloop;
  SeqDelay   trdelay;
  SeqObjList scan;
  SeqObjList dummypart;
  SeqObjList templatepart;
  SeqObjList grappapart;
  SeqObjList imagingpart;
  SeqObjList slicepart;
  SeqObjList slicepart_dummy;
  SeqObjList slicepart_template;
  SeqObjList slicepart_grappa;
  SeqObjList preppart;
  SeqObjList bladepart;
  SeqDelay   exc2acq;
  SeqDelay   exc2refoc;
  SeqDelay   refoc2acq;
  SeqGradConstPulse spoiler;
  SeqTrigger trigger;
  SeqGradTrapezParallel crusher;
  SeqDelay   crusherdelay;
  SeqFieldMap fmapscan;
  SeqVecIter phaseiter;

  SeqGradTrapez deph3d;
  SeqGradTrapez blip3d_read;
  SeqGradTrapez blip3d_phase;
  SeqGradTrapez blip3d_template_read;
  SeqGradTrapez blip3d_template_phase;
  SeqGradTrapez blip3d_grappa_read;
  SeqGradTrapez blip3d_grappa_phase;
  SeqGradVectorPulse blip3d_slice;
  SeqGradChanParallel blip3d;
  SeqGradChanParallel blip3d_template;
  SeqGradChanParallel blip3d_grappa;


  SeqObjLoop loop3d;
  SeqVector index3d;
  SeqVector epiindex;


  SeqObjList epiacq_part;
  SeqObjList epiacq_dummy_part;
  SeqObjList epiacq_template_part;
  SeqObjList epiacq_grappa_part;


  SeqRotMatrixVector bladerot;

  SeqMagnReset reset;


  dvector bladeangels;

public:

// This constructor creates an empty EPI sequence
METHOD_CLASS(const STD_string& label) : SeqMethod(label) {
  set_description("EPI-based Periodically Rotated Overlapping ParallEL Lines with Enhanced Reconstruction (PROPELLER) sequence. See for example J. G. Pipe, Magn. Reson. Med. 42:963-969 (1999)");
}

unsigned int numof_testcases() const {return 2;}


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);


  Blades=8;
  Blades.set_description("Number of PROPELLER blades");

  ShortAxis=false;
  ShortAxis.set_description("Readout direction along short axis of blades");

  BladeOversampling=0.0;
  BladeOversampling.set_minmaxval(0.0,100.0);
  BladeOversampling.set_unit("%").set_description("Oversampling (overlap) in short direction of blade");

  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");

  T1Ernst=0.0;
  T1Ernst.set_minmaxval(0.0,5000.0).set_description("If non-zero, the flip angle will be set to the Ernst angle using this T1 for optimum SNR");

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

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

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

  InvertPartialFourier=false;
  InvertPartialFourier.set_description("Invert position in phase encoding direction at which partial Fourier cut off is performed (i.e. beginning or end of readout)");

  UpDownBlade=false;
  UpDownBlade.set_description("Sample each blade twice with opposite phase-encoding direction");

  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");


  // register method parameters for user interface, parameter files, etc.
  append_parameter(Blades,"Blades");
  append_parameter(ShortAxis,"ShortAxis");

  append_parameter(EPIMode,"EPIMode");
  append_parameter(TakeMinEchoTime,"TakeMinEchoTime");
  append_parameter(DummyCycles,"DummyCycles");
  append_parameter(TemplateScan,"TemplateScan");
  append_parameter(RampSampling,"RampSampling");
  append_parameter(FatSaturation,"FatSaturation");
  append_parameter(T1Ernst,"T1Ernst");
  append_parameter(FieldMap,"FieldMap");

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

  append_parameter(UpDownBlade,"UpDownBlade");

  append_parameter(PulseDur,"PulseDur");

  append_parameter(fMRITrigger,"fMRITrigger");

  append_parameter(InvertPartialFourier,"InvertPartialFourier");

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


  // alternative default settings for sequence test
  if(get_current_testcase()==1) {
    geometryInfo->set_Mode(voxel_3d);
    ShortAxis=false;
    BladeOversampling=0.0;
    TemplateScan.set_actual("NoCorrection");
  }


}



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

  float gamma=systemInfo->get_gamma();

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


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

  // 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);

  if(geometryInfo->get_Mode()==voxel_3d) {
    slicethick=0.8*geometryInfo->get_FOV(sliceDirection);
    extra_slicethick_refoc=0.2*slicethick;
  }

  float spatres=slicethick/4.0;


  // calculate Ernst angle accordng to TR
  float flipangle=commonPars->get_FlipAngle();
  if(T1Ernst>0.0) {
    flipangle=180.0/PII * acos( exp ( -secureDivision ( commonPars->get_RepetitionTime(), T1Ernst) ) );
    commonPars->set_FlipAngle( flipangle );
  }

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


  // refocusing pulse
  refoc=SeqPulsarSinc("refoc",slicethick+extra_slicethick_refoc,false,PulseDur,180.0,spatres);
  refoc.set_freqlist( 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: //////////////////////////////

  // square FOV
  int sizeRadial=commonPars->get_MatrixSize(readDirection);
  commonPars->set_MatrixSize(phaseDirection,sizeRadial,noedit);

  int size3d=1;
  if(geometryInfo->get_Mode()==voxel_3d) {
    float resolution=secureDivision(geometryInfo->get_FOV(readDirection),commonPars->get_MatrixSize(readDirection));
    size3d=int(secureDivision(geometryInfo->get_FOV(sliceDirection),resolution)+0.5); // isotropic resolution
  }
  commonPars->set_MatrixSize(sliceDirection,size3d,noedit);


  int readpts_blade=sizeRadial;
  int pelines_blade=sizeRadial;

  int bladewidth=sizeRadial;

  if(Blades>1) {
    float min_bladewidth=sizeRadial*tan(0.5*PII/Blades);  // exact solution
    bladewidth=int( ( (1.0-0.01*BladeOversampling) * min_bladewidth + 0.01*BladeOversampling*sizeRadial ) +0.5 );
  }

  if(ShortAxis) readpts_blade=bladewidth;
  else          pelines_blade=bladewidth;


  float os_read=2.0; // For reduced undersampling artifacts

  float fov=geometryInfo->get_FOV(readDirection); // uniform FOV

  int epi_reduction=commonPars->get_ReductionFactor();
  float epi_partfour=commonPars->get_PartialFourier();

  if(geometryInfo->get_Mode()==voxel_3d) {
    epi_reduction=1;  // 3D phase encoding gets the reduction
    epi_partfour=0.0; // 3D phase encoding gets the partial fourier
  }


  epiacq=SeqAcqEPI("epiacq",commonPars->get_AcqSweepWidth(),
                   readpts_blade, fov,
                   pelines_blade, fov,
                   1, epi_reduction, os_read, "", 0, 0, rampType(int(RampMode)),
                   RampSampling, RampSteepness, epi_partfour, 0, InvertPartialFourier);

  // display sampling extents in read/phase direction
  NumOfGradEchoes=epiacq.get_numof_gradechoes();
  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);


  //////////////// 3D encoding: //////////////////////////

  if(geometryInfo->get_Mode()==voxel_3d) {


    // use SeqGradPhaseEnc as a template to calculate phase encoding
    SeqGradPhaseEnc pe3d("pe3d", commonPars->get_MatrixSize(sliceDirection), geometryInfo->get_FOV(sliceDirection), deph.get_gradduration(), sliceDirection,
                                 linearEncoding, noReorder, 1, commonPars->get_ReductionFactor(), DEFAULT_ACL_BANDS, commonPars->get_PartialFourier());


    fvector gradints3d=pe3d.get_trims()*pe3d.get_constduration()*pe3d.get_strength();
    ODINLOG(odinlog,normalDebug) << "gradints3d=" << gradints3d << STD_endl;

    const SeqVector& pe3dvec=pe3d; // access to vectorgrad
    ivector indexvec3d=pe3dvec.get_indexvec();
    ODINLOG(odinlog,normalDebug) << "indexvec3d=" << indexvec3d << STD_endl;


    float strength3d=0.7*systemInfo->get_max_grad();

    // dephaser before readout
    float dephintegral=0.0;
    if(gradints3d.size()) dephintegral=gradints3d[0];
    deph3d=SeqGradTrapez("deph3d", dephintegral, strength3d, sliceDirection);
    deph/=deph3d;
    deph_template/=deph3d;
    deph_grappa/=deph3d;

    // rephasers inbetween readouts
    fvector rewind_integral=-epiacq.get_gradintegral();
    fvector rewind_template_integral=-epiacq_template.get_gradintegral();
    fvector rewind_grappa_integral=-epiacq_grappa.get_gradintegral();
    blip3d_read= SeqGradTrapez("blip3d_read",  rewind_integral[readDirection],  strength3d, readDirection);
    blip3d_phase=SeqGradTrapez("blip3d_phase", rewind_integral[phaseDirection], strength3d, phaseDirection);
    blip3d_template_read= SeqGradTrapez("blip3d_template_read",  rewind_template_integral[readDirection],  strength3d, readDirection);
    blip3d_template_phase=SeqGradTrapez("blip3d_template_phase", rewind_template_integral[phaseDirection], strength3d, phaseDirection);
    blip3d_grappa_read= SeqGradTrapez("blip3d_grappa_read",  rewind_grappa_integral[readDirection],  strength3d, readDirection);
    blip3d_grappa_phase=SeqGradTrapez("blip3d_grappa_phase", rewind_grappa_integral[phaseDirection],  strength3d, phaseDirection);

    // blip in slice direction
    fvector blipintegrals(gradints3d.size()); blipintegrals=0.0;
    for(unsigned int i=0; i<(gradints3d.size()-1); i++) blipintegrals[i]=gradints3d[i+1]-gradints3d[i]; // differences for blips
    ODINLOG(odinlog,normalDebug) << "blipintegrals=" << blipintegrals << STD_endl;
    float maxabsintegral=blipintegrals.maxabs();
    float blipstrength_slice=sqrt(0.5*systemInfo->get_max_slew_rate()*maxabsintegral); // At least one half of the gradient as flat top
    float gradduration_slice=secureDivision(maxabsintegral, blipstrength_slice);
    blip3d_slice=SeqGradVectorPulse("blip3d_slice", sliceDirection, blipstrength_slice, blipintegrals/maxabsintegral, gradduration_slice);

/*
    if(UpDownBlade) {
      deph3d.invert_strength();
      blip3d_slice.invert_strength();
    }
*/

    blip3d=SeqGradChanParallel("blip3d");
    blip3d_template=SeqGradChanParallel("blip3d_template");
    blip3d_grappa=SeqGradChanParallel("blip3d_grappa");

    blip3d +=          blip3d_read          / blip3d_phase          / blip3d_slice;
    blip3d_template += blip3d_template_read / blip3d_template_phase / blip3d_slice;
    blip3d_grappa +=   blip3d_grappa_read   / blip3d_grappa_phase   / blip3d_slice;

    index3d=SeqVector("index3d");
    index3d.set_indexvec(indexvec3d);

    epiindex=SeqVector("epiindex",indexvec3d.size()); // to enumerate EPI readouts
  }


  /////////////////// Rotation of Blades ////////////////////////////////////////////////

  bladerot=SeqRotMatrixVector("bladerot");

  bool full_cycle=UpDownBlade; //(geometryInfo->get_Mode()!=voxel_3d && UpDownBlade);

  int nangles=Blades;
  if(full_cycle) nangles=2*Blades;

  bladeangels.resize(nangles);
  for(int iangle=0; iangle<nangles; iangle++) {
    RotMatrix rm("rotmatrix"+itos(iangle));
    int iblade=iangle;
    if(full_cycle) iblade=iangle/2;
    float ang=float(iblade)/float(Blades)*PII; // 0...180 deg
    if(full_cycle && iangle%2) ang+=PII; // adjacent blades have opposite direction
    rm.set_inplane_rotation(ang);
    bladeangels[iangle]=ang;
    bladerot.append(rm);
  }


  /////////////////// 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 blades
  bladeloop=SeqObjLoop("bladeloop");

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

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

  // loop to iterate over 3rd dimension phase encoding
  loop3d=SeqObjLoop("loop3d");

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

  trdelay=SeqDelay("trdelay");


  //////////////// 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",sliceDirection,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


  //////////////// trigger: //////////////////////////////

  trigger=SeqTrigger("fmri_trigger",1.0);

  reset=SeqMagnReset("reset");

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

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


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

  preppart+=reset;

  epiacq_part=SeqObjList("epiacq_part");
  epiacq_dummy_part=SeqObjList("epiacq_dummy_part");
  epiacq_template_part=SeqObjList("epiacq_template_part");
  epiacq_grappa_part=SeqObjList("epiacq_grappa_part");

  if(geometryInfo->get_Mode()==voxel_3d) {

    epiacq_part += loop3d( epiacq + blip3d )[index3d][epiindex][blip3d_slice];
    epiacq_dummy_part += loop3d( epiacq_dummy + blip3d )[index3d][epiindex][blip3d_slice];
    epiacq_template_part += loop3d( epiacq_template + blip3d_template )[index3d][epiindex][blip3d_slice];
    epiacq_grappa_part += loop3d( epiacq_grappa + blip3d_grappa )[index3d][epiindex][blip3d_slice];

  } else {

    epiacq_part += epiacq;
    epiacq_dummy_part += epiacq_dummy;
    epiacq_template_part += epiacq_template;
    epiacq_grappa_part += epiacq_grappa;

  }






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

  if(fMRITrigger && (geometryInfo->get_Mode()!=voxel_3d) ) slicepart+= trigger; // trigger for each single blade

  if (EPIMode=="FIDMode") {

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

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

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

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

  }

  if (EPIMode=="SEMode") {

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

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

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

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

  templatepart.set_gradrotmatrixvector(bladerot);
  grappapart.  set_gradrotmatrixvector(bladerot);
  imagingpart. set_gradrotmatrixvector(bladerot);


  if (EPIMode=="FIDMode") {
    slicepart_dummy =    sliceloop( dummypart    + trdelay )[exc];
    slicepart_template = sliceloop( templatepart + trdelay )[exc];
    slicepart_grappa =   sliceloop( grappapart   + trdelay )[exc];
    slicepart +=         sliceloop( imagingpart  + trdelay )[exc];
  }
  if (EPIMode=="SEMode") {
    slicepart_dummy =    sliceloop( dummypart + trdelay    )[exc][refoc];
    slicepart_template = sliceloop( templatepart + trdelay )[exc][refoc];
    slicepart_grappa  =  sliceloop( grappapart + trdelay   )[exc][refoc];
    slicepart +=         sliceloop( imagingpart + trdelay  )[exc][refoc];
  }


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

  if(FieldMap) scan += fmapscan + trdelay;


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


  if(TemplateScan=="PhaseCorrection") {
    scan += bladeloop(
              slicepart_template
            )[bladerot];
  }


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

  if(fMRITrigger && (geometryInfo->get_Mode()==voxel_3d) ) bladepart+= trigger; // trigger only for each full blade cycle

  bladepart+= bladeloop(
                slicepart
              )[bladerot];

  scan+= reploop( 
           bladepart
         )[commonPars->get_NumOfRepetitions()];


  set_sequence( scan );
}




void method_rels() {


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


   double acq_center=epiacq.get_acquisition_center();
   if(geometryInfo->get_Mode()==voxel_3d) acq_center= (1.0-commonPars->get_PartialFourier())*0.5*epiacq_part.get_duration();


   double min_echo_time=0.0;


   if (EPIMode=="FIDMode") {

     min_echo_time=(exc.get_duration()-exc.get_magnetic_center())+deph.get_duration()+acq_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()
                + acq_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
   trdelay=0.0;

   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
  epiacq.set_reco_vector(slice,exc);
  epiacq_template.set_reco_vector(slice,exc);
  epiacq_grappa.set_reco_vector(slice,exc);


  if(geometryInfo->get_Mode()==voxel_3d) {
    epiacq.set_reco_vector(line3d,index3d);
    epiacq_template.set_reco_vector(line3d,index3d);
    epiacq_grappa.set_reco_vector(line3d,index3d);

    epiacq.set_reco_vector(epi, epiindex);
    epiacq_template.set_reco_vector(epi, epiindex);
    epiacq_grappa.set_reco_vector(epi, epiindex);

    // For multi-frequency distortion correction
    int nepi=epiindex.get_vectorsize();
    dvector offsettime(nepi);
    double epidur=epiacq.get_duration() + blip3d.get_duration();
    for(int i=0; i<nepi; i++) offsettime[i]=i*epidur;
    recoInfo->set_DimValues(epi,offsettime);
  }


  epiacq.set_reco_vector(cycle,bladerot);
  epiacq_template.set_reco_vector(cycle,bladerot);
  epiacq_grappa.set_reco_vector(cycle,bladerot);
  recoInfo->set_DimValues(cycle,bladeangels);

//  if(geometryInfo->get_Mode()==voxel_3d && UpDownBlade) recoInfo->set_CmdLineOpts("-ip -sflip");

}

};


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


// entry point for the sequence module
ODINMETHOD_ENTRY_POINT