//////////////////////////////////////////////////////////////////////////////
//
// A simple spiral sequence
//
//////////////////////////////////////////////////////////////////////////////

// includes for the sequence objects
#include <odinseq/seqall.h>


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


// The method itself is a class that is derived from
// an abstract base class 'SeqMethod' that contains all
// routines common to all methods:

class METHOD_CLASS : public SeqMethod {

 public:

  // Constructor that takes the methods identifier (unique string) as its argument
  METHOD_CLASS(const STD_string& label);


  // virtual functions that are overwritten in this method to build the sequence,
  // calculate parameter relations, etc.:
  void method_pars_init();
  void method_seq_init();
  void method_rels();
  void method_pars_set();

 private:

  // Parameters for this method:
  LDRint        SpiralSegments;
  LDRbool       InOutSpiral;
  LDRfloat      OptimizedPar;

  LDRtrajectory Trajectory;
  LDRbool       FieldMap;
  LDRbool       FatSaturation;


  LDRfloatArr   DensityComp;
  LDRdouble     ReadoutDuration;
  LDRint        DisplaySegment;
  LDRfloatArr   kx;
  LDRfloatArr   ky;



  // Sequence objects for this method:
  SeqPulsar exc;
  SeqSat    fatsat;

  SeqDelay echodelay;
  SeqDelay relaxdelay;

  SeqObjLoop sliceloop;
  SeqObjLoop segloop;
  SeqObjLoop reploop;

  SeqObjList kernel;
  SeqObjList segpart;
  SeqObjList scan;

  SeqAcqSpiral spiral;

  SeqFieldMap fmapscan;

};


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

METHOD_CLASS::METHOD_CLASS (const STD_string& label)
                                   : SeqMethod(label) {
  // Put in here all stuff that will once be initialised and
  // never changed

  // Specify a short description of the method
  set_description("A simple sequence with different spiral k-space trajectories. "
                  "Optionally, a fieldmap can be acquired prior to the actual scan "
                  "which will be used for conjugate-phase multi-frequency reconstruction. ");


}



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

void METHOD_CLASS::method_pars_init() {
  // In this function the methods parameters will be initialised


  // Assign default values:
  commonPars->set_RepetitionTime(1000.0);
  commonPars->set_EchoTime(0.0);
  commonPars->set_AcqSweepWidth(100.0);
  commonPars->set_MatrixSize(readDirection,64);
  commonPars->set_MatrixSize(phaseDirection,64,noedit);

  Trajectory.set_function_mode(twoDeeMode);
  Trajectory.set_funcpars("WrapSpiral(32,0.80)"); // Default trajectory
  Trajectory.set_description("Spiral trajectory");
  append_parameter(Trajectory,"Trajectory");


  // Make method specific parameters visible in the user interface
  SpiralSegments=8;
  SpiralSegments.set_description("Number of spiral interleaves");
  append_parameter(SpiralSegments,"SpiralSegments");

  InOutSpiral=false;
  InOutSpiral.set_description("Spiral-in spiral-out readout");
  append_parameter(InOutSpiral,"InOutSpiral");

  OptimizedPar=-1.0;
  OptimizedPar.set_description("Optimize trajectory for minimum readout length if trajectory has a free tunable parameter. Set this parameter to a negative value to trigger optimization.");
  append_parameter(OptimizedPar,"OptimizedPar");

  ReadoutDuration.set_description("Total length of spiral readout");
  append_parameter(ReadoutDuration,"ReadoutDuration",noedit);


  DensityComp.set_filemode(exclude);
  DisplaySegment.set_filemode(exclude);
  kx.set_filemode(exclude);
  ky.set_filemode(exclude);


  if(systemInfo->get_platform()!=numaris_4) {
    append_parameter(DensityComp,"DensityComp");
    append_parameter(DisplaySegment,"DisplaySegment");
    append_parameter(kx,"kx");
    append_parameter(ky,"ky");
  }

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

  FieldMap=true;
  FieldMap.set_description("Fieldmap pre-scan for distortion correction");
  append_parameter(FieldMap,"FieldMap");

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


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


void METHOD_CLASS::method_seq_init() {
  // Put in here all stuff to create the layout of the sequence


  ///////////////// Excitation Pulse: /////////////////////

  // Get the slice thickness from the global geometry handler 'geometryInfo':
  float slicethick=geometryInfo->get_sliceThickness();


  // Create the sinc shaped excitation pulse.
  exc=SeqPulsarSinc("exc",slicethick,true,2.0,commonPars->get_FlipAngle());

  // Set the frequency list of the excitation pulse so that we will excite all slices in the slicepack
  exc.set_freqlist( systemInfo->get_gamma("") * exc.get_strength() / (2.0*PII) *  geometryInfo->get_sliceOffsetVector() );

  // This is useful for simulating/visualization of the sequence
  exc.set_pulse_type(excitation);


  // This loop object is used to loop over the slices
  sliceloop=SeqObjLoop("sliceloop");


  ///////////////// Fat saturation Pulse: /////////////////////

  fatsat=SeqSat("fatsat",fat);


  ////////////////// Geometry: /////////////////////////////////

  float fov=geometryInfo->get_FOV(readDirection);
  unsigned int sizeRadial=commonPars->get_MatrixSize(readDirection);
  commonPars->set_MatrixSize(phaseDirection,sizeRadial,noedit);

  // check bounds
  if(SpiralSegments<1) SpiralSegments=1;
  if(SpiralSegments>int(sizeRadial/2)) SpiralSegments=sizeRadial/2;


  //////////////// Spiral Readout: //////////////////////////////

  // Perform optimization only if requested manually, i.e. if OptimizedPar<0.0
  bool optimize=(OptimizedPar<0.0);

  if(!optimize) Trajectory.set_parameter(_TRAJ_OPTIMIZE_PARLABEL_,ftos(OptimizedPar));

  spiral=SeqAcqSpiral("spiral",commonPars->get_AcqSweepWidth(),fov,sizeRadial,SpiralSegments,Trajectory,InOutSpiral,optimize);

  // Cache calculated value (e.g. for calculation on VxWorks)
  if(optimize) OptimizedPar=atof(Trajectory.get_parameter(_TRAJ_OPTIMIZE_PARLABEL_).c_str());

  //////////////// Delays: //////////////////////////////

  // relaxation delay after each readout
  relaxdelay=SeqDelay("relaxdelay");

  // delay to obtain correct TE
  echodelay=SeqDelay("echodelay");


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

  // Construct a loop object to iterate through the phase encoding DisplaySegments
  segloop=SeqObjLoop("segloop");

  // Construct a loop object to perform repetitions of the experiment
  reploop=SeqObjLoop("reploop");


  //////////////// 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: ///////////////

  // This sequence container will hold the objects for one readout:
  kernel=SeqObjList("kernel");

  // This sequence container will hold the objects for one slicepack:
  segpart=SeqObjList("segpart");

  scan=SeqObjList("scan");


  if(FatSaturation) {
    kernel += fatsat;
  }


  kernel +=  exc  + echodelay  + spiral;

  segpart = segloop(
              sliceloop(
                kernel + relaxdelay
              )[exc]
            )[spiral.get_segment_vector()];


  if(FieldMap) scan += fmapscan + relaxdelay;

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


 // Finally, build the whole sequence
  set_sequence( scan );

}


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


void METHOD_CLASS::method_rels() {
  // Put in here all stuff that has to be performed whenever one of the sequence parameters
  // has been changed by the user


  // calculate relaxdelay to get the desired repetition time
  float scandur=kernel.get_duration()*float(geometryInfo->get_nSlices());
  if(scandur>commonPars->get_RepetitionTime()) commonPars->set_RepetitionTime(scandur);
  relaxdelay.set_duration( (commonPars->get_RepetitionTime()-scandur)/float(geometryInfo->get_nSlices()) );

  // set TE
  float min_echo_time=(exc.get_duration()-exc.get_magnetic_center())+(spiral.get_acquisition_center());
  if (commonPars->get_EchoTime()<min_echo_time) commonPars->set_EchoTime(min_echo_time);
  echodelay=commonPars->get_EchoTime()-min_echo_time;



  ReadoutDuration=secureDivision(double(spiral.get_npts()),commonPars->get_AcqSweepWidth()) ;

  if(systemInfo->get_platform()!=numaris_4) {
    DensityComp=spiral.get_denscomp();
    kx=spiral.get_ktraj(DisplaySegment,readDirection);
    ky=spiral.get_ktraj(DisplaySegment,phaseDirection);
  }


}


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

void METHOD_CLASS::method_pars_set() {
  Log<Seq> odinlog(this,"method_pars_set");

  // Put in here all stuff that has to be performed after the parameters have been edited by the user
  // and before the sequence is played out

  spiral.set_reco_vector(slice,exc);
}


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

// entry point for the sequence module
ODINMETHOD_ENTRY_POINT