C++ dbl函数代码示例

pulsesequence()
{
   char         Cdseq[MAXSTR], refoc[MAXSTR], sspul[MAXSTR];
   int          mult = (0.5 + getval("mult"));
   double	rg1 = 2.0e-6,
                j1xh = getval("j1xh"),
         	gt0 = getval("gt0"),
         	gzlvl0 = getval("gzlvl0"),
         	pp = getval("pp"),
	 	pplvl = getval("pplvl"),
	 	compH = getval("compH"),
		Cdpwr = getval("Cdpwr"),
	 	Cdres = getval("Cdres"), 
	 	tau1 = 0.002,
	 	tau2 = 0.002,
	 	Cdmf = getval("Cdmf");
   
   shape   hdx;

/* Get new variables from parameter table */

   getstr("Cdseq", Cdseq);
   getstr("refoc", refoc);
   getstr("sspul", sspul);
   
   if (j1xh < 1.0) j1xh = 150.0;

   hdx = pboxHT_F1i("gaus180", pp*compH, pplvl); /* HADAMARD stuff */

   if (getval("htcal1") > 0.5)          /* Optional fine power calibration */
     hdx.pwr = getval("htpwr1");

   if(j1xh > 0.0)
     tau1 = 0.25/j1xh;
   tau2 = tau1;
   if (mult > 2) tau2=tau1/2.5;
   else if ((mult == 0) || (mult == 2)) tau2=tau1/2.0;

   dbl(ct, v1);                     /*  v1 = 02 */
   add(two,v1,oph);
   mod4(oph,oph);                   /* oph = 02 */


/* Calculate delays */

   if (dm[0] == 'y')
   {
     abort_message("decoupler must be set as dm=nny\n");
   }
   if(refoc[A]=='n' && dm[C]=='y')
   {
     abort_message("with refoc=n decoupler must be set as dm=nnn\n");
   }
 
/* Relaxation delay */

   status(A);

   delay(0.05);
   zgradpulse(gzlvl0,gt0);
   if (sspul[0] == 'y')
   {
     obspower(tpwr); decpower(pplvl);
     simpulse(pw, pp, zero, zero, rg1, rg1);  /* destroy H and C magnetization */
     zgradpulse(gzlvl0,gt0);
   }

   delay(d1);
 
   status(B);

     obspower(Cdpwr);
     obsunblank(); xmtron(); 
     obsprgon(Cdseq, 1.0/Cdmf, Cdres);  

     pbox_decpulse(&hdx, zero, rg1, rg1);

     obsprgoff(); xmtroff(); obsblank();  
     obspower(tpwr); decpower(pplvl);
          
     delay(2.0e-4);		
     
     decrgpulse(pp, v1, rg1, rg1);
     
     delay(tau1 - POWER_DELAY);
     simpulse(2.0*pw, 2.0*pp, zero, zero, rg1, rg1);
     txphase(one); decphase(one); 
     delay(tau1);
     simpulse(pw, pp, one, one, rg1, rg1);
     if(refoc[A]=='y')
     {
       txphase(zero); decphase(zero); 
       delay(tau2);
       simpulse(2.0*pw, 2.0*pp, zero, zero, rg1, rg1);
       delay(tau2 - rof2 - POWER_DELAY);
       decrgpulse(pp, zero, rg1, rof2);
     }
     decpower(dpwr);    

   status(C);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  /* Calculate total acquisition time */
  g_setExpTime(tr*(nt*nseg*getval("arraydim") + ssc) + trimage*getval("arraydim"));


  /* Return parameters to VnmrJ */
  putvalue("rgss",ss_grad.tramp);  //90  slice ramp
  if (ss2_grad.enableButterfly) {   //180 slice ramps
    putvalue("rcrush",ss2_grad.crusher1RampToCrusherDuration);
    putvalue("rgss2",ss2_grad.crusher1RampToSsDuration);
  }
  else {
    putvalue("rgss2",ss2_grad.tramp);
  }
  if (ro_grad.enableButterfly) {
    putvalue("rgro",ro_grad.crusher1RampToSsDuration);
  }
  else {   
    putvalue("rgro",ro_grad.tramp);      //RO ramp
  }
  putvalue("tror",ror_grad.duration);  //ROR duration
  putvalue("rgror",ror_grad.tramp);    //ROR ramp
  putvalue("gpe",pe_grad.peamp);         //PE max amp
  putvalue("gss",ss_grad.ssamp);
  putvalue("gro",ro_grad.roamp);



  /* PULSE SEQUENCE *************************************/
  initval(fabs(ssc),vssc);      // Compressed steady-state counter
  assign(one,vacquire);         // real-time acquire flag

  /* Phase cycle: Alternate 180 phase to cancel residual FID */
  mod2(ct,vphase180);           // 0101
  dbl(vphase180,vphase180);     // 0202
  add(vphase180,one,vphase180); // 1313 Phase difference from 90
  add(vphase180,oph,vphase180);

  obsoffset(resto);
  delay(4e-6);
    
  initval(nseg,vseg);
  loop(vseg,vseg_ctr);

    /* TTL scope trigger **********************************/       
    sp1on(); delay(4e-6); sp1off();

    /* Compressed steady-states: 1st array & transient, all arrays if ssc is negative */
    if ((ix > 1) && (ssc > 0))
      assign(zero,vssc);
    sub(vseg_ctr,vssc,vseg_ctr);   // vpe_ctr counts up from -ssc
    assign(zero,vssc);
    ifzero(vseg_ctr);
      assign(zero,vacquire);       // Start acquiring when vseg_ctr reaches zero
    endif(vseg_ctr);
    setacqvar(vacquire);           // Turn on acquire when vacquire is zero

    if (ticks) {
      xgate(ticks);
      grad_advance(gpropdelay);
      delay(4e-6);
    }

    if(ir[0] == 'y') {  /* IR for all slices prior to data acquisition */
      obspower(ir_rf.powerCoarse);
      obspwrf(ir_rf.powerFine);
      delay(4e-6);

void pulsesequence()
{
   double	   slpwrR = getval("slpwrR"),
		   slpwR = getval("slpwR"),
		   mixR = getval("mixR"),
                   selpwrA = getval("selpwrA"),
                   selpwA = getval("selpwA"),
                   gzlvlA = getval("gzlvlA"),
                   gtA = getval("gtA"),
                   selpwrB = getval("selpwrB"),
                   selpwB = getval("selpwB"),
                   gzlvlB = getval("gzlvlB"),
                   gtB = getval("gtB"),
                   gstab = getval("gstab"),
		   selfrq = getval("selfrq"),
                   gzlvl1 = getval("gzlvl1"),
                   gt1 = getval("gt1"),
                   gzlvl2 = getval("gzlvl2"),
                   gt2 = getval("gt2"),
                   zqfpw1 = getval("zqfpw1"),
                   zqfpwr1 = getval("zqfpwr1"),
                   gzlvlzq1 = getval("gzlvlzq1"),
                   phincr1 = getval("phincr1");
   char            selshapeA[MAXSTR],selshapeB[MAXSTR], slpatR[MAXSTR],
                   zqfpat1[MAXSTR], alt_grd[MAXSTR];

//synchronize gradients to srate for probetype='nano'
//   Preserve gradient "area"
        gtA = syncGradTime("gtA","gzlvlA",1.0);
        gzlvlA = syncGradLvl("gtA","gzlvlA",1.0);
        gtB = syncGradTime("gtB","gzlvlB",1.0);
        gzlvlB = syncGradLvl("gtB","gzlvlB",1.0);

   getstr("slpatR",slpatR);
   getstr("selshapeA",selshapeA);
   getstr("selshapeB",selshapeB);
   getstr("zqfpat1",zqfpat1);
   getstr("alt_grd",alt_grd);

   if (strcmp(slpatR,"cw") &&
        strcmp(slpatR,"troesy") &&
        strcmp(slpatR,"dante"))
        abort_message("SpinLock pattern %s not supported!.\n", slpatR);

/* STEADY-STATE PHASECYCLING */
/* This section determines if the phase calculations trigger off of (SS - SSCTR) or off of CT */

  assign(ct,v17);
   ifzero(ssctr);
      assign(v17,v13);
   elsenz(ssctr);
                /* purge option does not adjust v13 during steady state */
      sub(ssval, ssctr, v13);
   endif(ssctr);

/* Beginning phase cycling */

   dbl(v13,v1);		/* v1 = 0 2 */
   hlv(v13,v13);
   dbl(v13,v20);		/* v20 = 00 22 */
   hlv(v13,v13);
   dbl(v13,v6);		/* v6 = 0000 2222 */
   hlv(v13,v13);
   dbl(v13,v7);		/* v7 = 00000000 22222222 */

   assign(v1,oph);

   if (getflag("Gzqfilt"))
      add(v7,oph,oph);

/* CYCLOPS */

   assign(v13,v14);	/* v14 = 8x0 8x1 8x2 8x3 */
   
   if (getflag("Gzqfilt"))
      hlv(v13,v14);	/* v14 = 16x0 16x1 16x2 16x3 */

   add(v1,v14,v1);      
   add(v20,v14,v20);      
   add(v6,v14,v6);      
   add(v7,v14,v7);      
   add(oph,v14,oph);

/*  add(oph,v18,oph);
  add(oph,v19,oph); */
  assign(zero,v9);

   mod2(ct,v2);    /* 01 01 */
   hlv(ct,v11); hlv(v11,v11); mod2(v11,v11); dbl(v11,v11); /* 0000 2222 */
   add(v11,v2,v11); mod4(v11,v11); /* 0101 2323  first echo in Excitation Sculpting */
   hlv(ct,v4); mod2(v4,v4);    /* 0011 */
   hlv(ct,v12); hlv(v12,v12); hlv(v12,v12); dbl(v12,v12); add(v12,v4,v12);
   mod4(v12,v12);   /* 0011 0011 2233 2233 second echo in Excitation Sculpting */

   dbl(v2,v2);    /* 0202 */
   dbl(v4,v4);    /* 0022 */
   add(v2,v4,v4); /* 0220 correct oph for Excitation Sculpting */
   add(oph,v4,oph); mod4(oph,oph);

   if (!strcmp(slpatR,"troesy")) 
//.........这里部分代码省略.........

void pulsesequence()
{
   char    c1d[MAXSTR];               /* option to record only 1D C13 spectrum */
   int     ncyc;
   double  tau1 = 0.002,         			          /*  t1 delay */
           post_del = 0.0,
           pwClvl = getval("pwClvl"), 	  /* coarse power for C13 pulse */
           pwC = getval("pwC"),     	  /* C13 90 degree pulse length at pwClvl */
	   compC = getval("compC"),
	   compH = getval("compH"),
	   mixpwr = getval("mixpwr"),
	   jCH = getval("jCH"),
	   gt0 = getval("gt0"),  		       
	   gt1 = getval("gt1"),  		       
	   gt2 = getval("gt2"),  		       
	   gzlvl0 = getval("gzlvl0"),
	   gzlvl1 = getval("gzlvl1"),
	   gzlvl2 = getval("gzlvl2"),
	   grec = getval("grec"),
	   phase = getval("phase");

           getstr("c1d",c1d);
	   
           ncyc=1;
           if(jCH > 0.0)
             tau1 = 0.25/jCH;

           dbl(ct, v1);                     /* v1 = 0 */
           mod4(v1,oph);
           hlv(ct,v2);
           add(v2,v1,v2);
           if (phase > 1.5)
             incr(v1);                      /* hypercomplex phase increment */           

           initval(2.0*(double)((int)(d2*getval("sw1")+0.5)%2),v10); 
           add(v1,v10,v1);
           add(oph,v10,oph);
           mod4(v1,v1);  mod4(v2,v2); mod4(oph,oph); 
           assign(zero,v3);

           if(FIRST_FID) 
	   {
	     HHmix = pbox_mix("HHmix", "DIPSI2", mixpwr, pw*compH, tpwr);  
	     if(c1d[A] == 'n')
	     {
	       opx("CHdec"); setwave("WURST2 30k/1.2m"); pbox_par("steps","600"); cpx(pwC*compC, pwClvl);
	       CHdec = getDsh("CHdec");
	     }
	   }
	   ncyc = (int) (at/HHmix.pw) + 1;
	   post_del = ncyc*HHmix.pw - at;
	            
             
/* BEGIN PULSE SEQUENCE */

      status(A);

	zgradpulse(gzlvl0, gt0);
	rgpulse(pw, zero, 0.0, 0.0);  /* destroy H-1 magnetization*/
	zgradpulse(gzlvl0, gt0);
	delay(1.0e-4);
	obspower(tpwr);
	txphase(v1);
        decphase(zero);
        dec2phase(zero);

        presat();
	obspower(tpwr);
        	
	delay(1.0e-5);

      status(B);

        if(c1d[A] == 'y')
	{
   	  rgpulse(pw,v1,0.0,0.0);                 /* 1H pulse excitation */
          delay(d2);
   	  rgpulse(pw,two,0.0,0.0);                 /* 1H pulse excitation */
          assign(oph,v3);
	}
	else
	{
          decunblank(); pbox_decon(&CHdec);

   	  rgpulse(pw,v1,0.0,0.0);                 /* 1H pulse excitation */
   	  txphase(zero);
   	
          delay(d2);

          pbox_decoff(); decblank();  
          decpower(pwClvl); decpwrf(4095.0);
   	  
	  delay(tau1 - POWER_DELAY);
          simpulse(2.0*pw, 2.0*pwC, zero, zero, 0.0, 0.0);
          txphase(one); decphase(one); dec2phase(one);
	  delay(tau1);
          simpulse(pw, pwC, one, one, 0.0, 0.0);
          txphase(zero); decphase(zero); dec2phase(zero);
	  delay(tau1);
          simpulse(2.0*pw, 2.0*pwC, zero, zero, 0.0, 0.0);
//.........这里部分代码省略.........

pulsesequence()

{

   double   gzlvlE = getval("gzlvlE"),
            gtE = getval("gtE"),
            EDratio = getval("EDratio"),
            gstab = getval("gstab"),
            mult = getval("mult"),
	    pwx180 = getval("pwx180"),
	    pwxlvl180 = getval("pwxlvl180"),
	    pwx180r = getval("pwx180r"),
	    pwxlvl180r = getval("pwxlvl180r"),
            tpwr180 = getval("tpwr180"),
            pw180 = getval("pw180"),
            hsglvl = getval("hsglvl"),
            hsgt = getval("hsgt"),
            gzlvl5 = getval("gzlvl5"),
            tau,
            taug;
   int      icosel,
	    phase1 = (int)(getval("phase")+0.5),
            prgcycle = (int)(getval("prgcycle")+0.5),
	    ZZgsign;
   char	    pwx180ad[MAXSTR],
	    pwx180ref[MAXSTR],
	    bipflg[MAXSTR],
	    pw180ad[MAXSTR];

//synchronize gradients to srate for probetype='nano'
//   Preserve gradient "area"
        gtE = syncGradTime("gtE","gzlvlE",0.5);
        gzlvlE = syncGradLvl("gtE","gzlvlE",0.5);

  getstr("pwx180ad", pwx180ad);
  getstr("pwx180ref", pwx180ref);
  getstr("pw180ad", pw180ad);
  getstr("bipflg",bipflg);
  if (bipflg[0] == 'y')
  {
        tpwr180=getval("tnbippwr");
        pw180=getval("tnbippw");
        getstr("tnbipshp",pw180ad);
  }
  if (bipflg[1] == 'y')
  {
        pwxlvl180=getval("dnbippwr");
        pwx180=getval("dnbippw");
        getstr("dnbipshp",pwx180ad);
  }

  tau = 1 / (4*(getval("j1xh")));

  if (mult > 0.5)
    taug = 2*tau + getval("tauC");
  else
    taug = gtE + gstab + 2 * GRADIENT_DELAY;
  ZZgsign=-1; 
  if (mult == 2) ZZgsign=1;
  icosel = 1;
 
  assign(ct,v17);
  assign(zero,v18);
  assign(zero,v19);

  if (getflag("prgflg") && (satmode[0] == 'y') && (prgcycle > 1.5))
    {
        hlv(ct,v17);
        mod2(ct,v18); dbl(v18,v18);
        if (prgcycle > 2.5)
           {
                hlv(v17,v17);
                hlv(ct,v19); mod2(v19,v19); dbl(v19,v19);
           }
     }

  settable(t1, 4, ph1);
  settable(t2, 2, ph2);
  settable(t3, 8, ph3);
  settable(t4, 16, ph4);
  settable(t5, 16, ph5);

  getelem(t1, v17, v1);
  getelem(t3, v17, v3);
  getelem(t4, v17, v4);
  getelem(t2, v17, v2);
  getelem(t5, v17, oph);

  assign(zero,v6);
  add(oph,v18,oph);
  add(oph,v19,oph);

  if ((phase1 == 2) || (phase1 == 5))
    icosel = -1;

  initval(2.0*(double)(((int)(d2*getval("sw1")+0.5)%2)),v14);
/*
  mod2(id2, v14);
  dbl(v14,v14);
*/
//.........这里部分代码省略.........

		// Return value of constant
		dbl FreshEval_d(std::shared_ptr<Variable> const& diff_variable) const override
		{
			return dbl(double(true_value_real_),double(true_value_imag_));
		}

pulsesequence() 
{ 
 double gstab = getval("gstab"),
	gt1 = getval("gt1"),
	gzlvl1 = getval("gzlvl1"),
	gt2 = getval("gt2"),
	gzlvl2 = getval("gzlvl2"),
        satpwr = getval("satpwr"),
        satdly = getval("satdly"),
	del = getval("del"),
	del2 = getval("del2"),
	dosyfrq = getval("sfrq"),
        gzlvlhs = getval("gzlvlhs"),
        hsgt = getval("hsgt"),
	Ddelta,dosytimecubed; 
 char convcomp[MAXSTR],satmode[MAXSTR],alt_grd[MAXSTR],lkgate_flg[MAXSTR],
      sspul[MAXSTR]; 
 
 getstr("convcomp",convcomp); 
 getstr("satmode",satmode);
 getstr("alt_grd",alt_grd);
 getstr("lkgate_flg",lkgate_flg);
 getstr("sspul",sspul);

       //synchronize gradients to srate for probetype='nano'
       //   Preserve gradient "area"
          gt1 = syncGradTime("gt1","gzlvl1",1.0);
          gzlvl1 = syncGradLvl("gt1","gzlvl1",1.0);
          gt2 = syncGradTime("gt2","gzlvl2",1.0);
          gzlvl2 = syncGradLvl("gt2","gzlvl2",1.0);

/* CHECK CONDITIONS */ 
  if (p1 == 0.0) p1 = 2*pw;
 
/* STEADY-STATE PHASECYCLING 
 This section determines if the phase calculations trigger off of (SS - SSCTR) 
   or off of CT */ 
 
   ifzero(ssctr); 
      dbl(ct, v1); 
      hlv(ct, v3); 
   elsenz(ssctr); 
      sub(ssval, ssctr, v7);	/* v7 = 0,...,ss-1 */ 
      dbl(v7, v1); 
      hlv(v7, v3); 
   endif(ssctr); 
 
/* PHASECYCLE CALCULATION */ 
   hlv(v3, v2); 
   mod2(v3, v3); 
   add(v3, v1, v1); 
   assign(v1, oph); 
   dbl(v2, v4); 
   add(v4, oph, oph); 
   add(v2, v3, v2); 
 
 	Ddelta=gt1;             /*the diffusion-encoding pulse width is gt1*/ 
	if (convcomp[A]=='y') 
		dosytimecubed=Ddelta*Ddelta*(del - (4.0*Ddelta/3.0)); 
        else 
		dosytimecubed=Ddelta*Ddelta*(del - (Ddelta/3.0)); 
 
   putCmd("makedosyparams(%e,%e)\n",dosytimecubed,dosyfrq); 
   mod2(ct,v10);        /* gradients change sign at odd transients */
 
/* BEGIN ACTUAL SEQUENCE */ 
 status(A); 
     if (sspul[0]=='y')
       {
         zgradpulse(gzlvlhs,hsgt);
         rgpulse(pw,zero,rof1,rof1);
         zgradpulse(gzlvlhs,hsgt);
       }
     if (satmode[0] == 'y')
      {
        if (d1 - satdly > 0) delay(d1 - satdly);
        obspower(satpwr);
        rgpulse(satdly,zero,rof1,rof1);
        obspower(tpwr);
      }
     else delay(d1);

 if (getflag("wet")) wet4(zero,one);

 status(B); 
   if (del>0.0) { 
      if (convcomp[A]=='y')
	{ 
                if (lkgate_flg[0] == 'y')  lk_hold(); /* turn lock sampling off */
                if (alt_grd[0] == 'y')
                 { ifzero(v10);
                     zgradpulse(-gzlvl2,2.0*gt2);
                   elsenz(v10);
                     zgradpulse(gzlvl2,2.0*gt2);
                   endif(v10);
                 }
                else zgradpulse(-gzlvl2,2.0*gt2);
		delay(gstab); 
             rgpulse(pw, v1, rof1, rof1); 
		delay(d2/2.0); 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  }
  else {   
    putvalue("rgro",ro_grad.tramp);      //RO ramp
  }
  putvalue("tror",ror_grad.duration);  //ROR duration
  putvalue("rgror",ror_grad.tramp);    //ROR ramp
  putvalue("gpe",pe_grad.peamp);         //PE max amp
  putvalue("gss",ss_grad.ssamp);
  putvalue("gro",ro_grad.roamp);


  g_setExpTime(tr*(nt*pe_steps*arraydim + ssc));


  /* PULSE SEQUENCE *************************************/
  rotate();
  obsoffset(resto);
  roff = -poffset(pro,ro_grad.roamp);
  delay(4e-6);

  /* Begin phase-encode loop ****************************/       
  peloop(seqcon[2],pe_steps,vpe_steps,vpe_ctr);

    /* Read external kspace table if set ******************/       
    if (table)
      getelem(t1,vpe_ctr,vpe_index);
    else
      sub(vpe_ctr,vpe_offset,vpe_index);

    settable(t2,2,ph180);        // initialize phase tables and variables
    getelem(t2,vpe_ctr,vph180);
    add(oph,vph180,vph180);      // 180 deg pulse phase alternates +/- 90 from receiver
    mod2(ct,vph2);
    dbl(vph2,vph2);
    add(vph180,vph2,vph180);     // Alternate phase for 180 on odd transients

    /* Begin multislice loop ******************************/       
    msloop(seqcon[1],ns,vms_slices,vms_ctr);
      if (ticks) {
        xgate(ticks);
        grad_advance(gpropdelay);
      }
      /* TTL scope trigger **********************************/       
       sp1on(); delay(5e-6); sp1off();

      /* Prepulses ******************************************/       
      if (sat[0]  == 'y') satbands();
      if (fsat[0] == 'y') fatsat();
      if (mt[0]   == 'y') mtc();

      /* Optional IR pulse **********************************/ 
      if (ir[0] == 'y') {
	obspower(ir_rf.powerCoarse);
	obspwrf(ir_rf.powerFine);
	delay(4e-6);
	obl_shapedgradient(ssi_grad.name,ssi_grad.duration,0,0,ssi_grad.amp,NOWAIT);
	delay(ssi_grad.rfDelayFront);
	shapedpulselist(shapeIR,ssi_grad.rfDuration,oph,rof1,rof2,seqcon[1],vms_ctr);
	delay(ssi_grad.rfDelayBack);
	delay(ti_delay);
      }

      /* Slice select RF pulse ******************************/ 
      obspower(p1_rf.powerCoarse);
      obspwrf(p1_rf.powerFine);
      delay(4e-6);

	static inline void add(EcT& R, const EcT& P, const EcT& Q)
	{
		if (P.isZero()) { R = Q; return; }
		if (Q.isZero()) { R = P; return; }
#if MIE_EC_COORD == MIE_EC_USE_JACOBI
		Fp r, U1, S1, H, H3;
		Fp::square(r, P.z);
		Fp::square(S1, Q.z);
		Fp::mul(U1, P.x, S1);
		Fp::mul(H, Q.x, r);
		H -= U1;
		r *= P.z;
		S1 *= Q.z;
		S1 *= P.y;
		Fp::mul(r, Q.y, r);
		r -= S1;
		if (H.isZero()) {
			if (r.isZero()) {
				dbl(R, P, false);
			} else {
				R.clear();
			}
			return;
		}
		Fp::mul(R.z, P.z, Q.z);
		R.z *= H;
		Fp::square(H3, H); // H^2
		Fp::square(R.y, r); // r^2
		U1 *= H3; // U1 H^2
		H3 *= H; // H^3
		R.y -= U1;
		R.y -= U1;
		Fp::sub(R.x, R.y, H3);
		U1 -= R.x;
		U1 *= r;
		H3 *= S1;
		Fp::sub(R.y, U1, H3);
#elif MIE_EC_COORD == MIE_EC_USE_PROJ
		Fp r, PyQz, v, A, vv;
		Fp::mul(r, P.x, Q.z);
		Fp::mul(PyQz, P.y, Q.z);
		Fp::mul(A, Q.y, P.z);
		Fp::mul(v, Q.x, P.z);
		v -= r;
		if (v.isZero()) {
			Fp::add(vv, A, PyQz);
			if (vv.isZero()) {
				R.clear();
			} else {
				dbl(R, P, false);
			}
			return;
		}
		Fp::sub(R.y, A, PyQz);
		Fp::square(A, R.y);
		Fp::square(vv, v);
		r *= vv;
		vv *= v;
		Fp::mul(R.z, P.z, Q.z);
		A *= R.z;
		R.z *= vv;
		A -= vv;
		vv *= PyQz;
		A -= r;
		A -= r;
		Fp::mul(R.x, v, A);
		r -= A;
		R.y *= r;
		R.y -= vv;
#else
		Fp t;
		Fp::neg(t, Q.y);
		if (P.y == t) { R.clear(); return; }
		Fp::sub(t, Q.x, P.x);
		if (t.isZero()) {
			dbl(R, P, false);
			return;
		}
		Fp s;
		Fp::sub(s, Q.y, P.y);
		Fp::div(t, s, t);
		R.inf_ = false;
		Fp x3;
		Fp::square(x3, t);
		x3 -= P.x;
		x3 -= Q.x;
		Fp::sub(s, P.x, x3);
		s *= t;
		Fp::sub(R.y, s, P.y);
		R.x = x3;
#endif
	}

void pulsesequence()

{

   double   selpwrA = getval("selpwrA"),
            selpwA = getval("selpwA"),
            gzlvlA = getval("gzlvlA"),
            gtA = getval("gtA"),
            selpwrB = getval("selpwrB"),
            selpwB = getval("selpwB"),
            gzlvlB = getval("gzlvlB"),
            gtB = getval("gtB"),
   	    gzlvlE = getval("gzlvlE"),
            gtE = getval("gtE"),
            EDratio = getval("EDratio"),
            gstab = getval("gstab"),
	    pwx180 = getval("pwx180"),
	    pwxlvl180 = getval("pwxlvl180"),
            hsglvl = getval("hsglvl"),
            hsgt = getval("hsgt"),
            tau,
	    evolcorr,
	    jdly,
	    bird;
   int      icosel,
	    ijscale,
            prgcycle = (int)(getval("prgcycle")+0.5),
	    phase1=(int)(getval("phase")+0.5);
   char     selshapeA[MAXSTR],
            selshapeB[MAXSTR],
   	    pwx180ad[MAXSTR],
	    pwx180adR[MAXSTR];

//synchronize gradients to srate for probetype='nano'
//   Preserve gradient "area"
        gtE = syncGradTime("gtE","gzlvlE",0.5);
        gzlvlE = syncGradLvl("gtE","gzlvlE",0.5);
        gtA = syncGradTime("gtA","gzlvlA",1.0);
        gzlvlA = syncGradLvl("gtA","gzlvlA",1.0);
        gtB = syncGradTime("gtB","gzlvlB",1.0);
        gzlvlB = syncGradLvl("gtB","gzlvlB",1.0);

  getstr("selshapeA",selshapeA);
  getstr("selshapeB",selshapeB);
  getstr("pwx180ad", pwx180ad);
  getstr("pwx180adR", pwx180adR);
  evolcorr=(4*pwx/PI)+2*pw+8.0e-6;
  bird = 1 / (4*(getval("j1xh")));
  tau = 1/(4*(getval("jnxh")));
  ijscale=(int)(getval("jscale") + 0.5);
  jdly=(d2*ijscale);

  icosel = 1;

  assign(ct,v17);
  assign(zero,v18);
  assign(zero,v19);

  if (getflag("prgflg") && (satmode[0] == 'y') && (prgcycle > 1.5))
    {
        hlv(ct,v17);
        mod2(ct,v18); dbl(v18,v18);
        if (prgcycle > 2.5)
           {
                hlv(v17,v17);
                hlv(ct,v19); mod2(v19,v19); dbl(v19,v19);
           }
     }

   settable(t1,4,ph1);
   settable(t2,2,ph2);
   settable(t4,8,ph4);
   settable(t5,8,ph5);

   getelem(t1,v17,v1);
   getelem(t2,v17,v2);
   getelem(t4,v17,v4);
   getelem(t5,v17,oph);

  add(oph,v18,oph);
  add(oph,v19,oph);

  if ((phase1 == 2) || (phase1 == 5))
    icosel = -1;

  initval(2.0*(double)(((int)(d2*getval("sw1")+0.5)%2)),v14);
/*
  mod2(id2, v14);
  dbl(v14,v14);
*/
  add(v2, v14, v2);
  add(oph, v14, oph);

status(A);
   obspower(tpwr);

   delay(5.0e-5);
   if (getflag("sspul"))
        steadystate();
   if (satmode[0] == 'y')
//.........这里部分代码省略.........

void pulsesequence()
{
   double          p1lvl,
                   trim,
                   mix,
                   window,
                   cycles;
   char            sspul[MAXSTR];


/* LOAD AND INITIALIZE VARIABLES */
   mix = getval("mix");
   trim = getval("trim");
   p1lvl = getval("p1lvl");
   window=getval("window");
   getstr("sspul", sspul);

/* CHECK CONDITIONS */
   if ((phase1 != 3) && (arrayelements > 2))
   {
      fprintf(stdout, "PHASE=3 is required if MIX is arrayed!\n");
      psg_abort(1);
   }
   if (satdly > 9.999)
   {
      fprintf(stdout, "Presaturation period is too long.\n");
      psg_abort(1);
   }
   if (!newtransamp)
   {
      fprintf(stdout, "TOCSY requires linear amplifiers on transmitter.\n");
      fprintf(stdout, "Use DECTOCSY with the appropriate re-cabling,\n");
      psg_abort(1);
   }
   if ((p1 == 0.0) && (ix == 1))
      fprintf(stdout, "Warning:  P1 has a zero value.\n");
   if ((rof1 < 9.9e-6) && (ix == 1))
      fprintf(stdout,"Warning:  ROF1 is less than 10 us\n");

   if (satpwr > 40)
        {
         printf("satpwr too large  - acquisition aborted./n");
         psg_abort(1);
        }

/* STEADY-STATE PHASECYCLING */
/* This section determines if the phase calculations trigger off of (SS - SSCTR)
   or off of CT */

   ifzero(ssctr);
      hlv(ct, v13);
      mod2(ct, v1);
      hlv(ct, v2);
   elsenz(ssctr);
      sub(ssval, ssctr, v12);	/* v12 = 0,...,ss-1 */
      hlv(v12, v13);
      mod2(v12, v1);
      hlv(v12, v2);
   endif(ssctr);


/* CALCULATE PHASES */
/* A 2-step cycle is performed on the first pulse (90 degrees) to suppress
   axial peaks in the first instance.  Second, the 2-step F2 quadrature image
   suppression subcycle is added to all pulse phases and receiver phase.
   Finally, a 2-step cycle is performed on the spin-lock pulses. */

   mod2(v13, v13);
   dbl(v1, v1);
   incr(v1);
   hlv(v2, v2);
   mod2(v2, v2);
   dbl(v2, v2);
   incr(v2);
   add(v13, v2, v2);
   sub(v2, one, v3);
   add(two, v2, v4);
   add(two, v3, v5);
   add(v1, v13, v1);
   assign(v1, oph);
   if (phase1 == 2)
      incr(v1);
   if (phase1 == 3)
      add(v1, id2, v1);

 /*HYPERCOMPLEX MODE USES REDFIELD TRICK TO MOVE AXIAL PEAKS TO EDGE*/
   if ((phase1==1)||(phase1==2))
   {
      initval(2.0*(double)(d2_index%2),v6);
      add(v1,v6,v1); add(oph,v6,oph);
   } 

/* CALCULATE AND INITIALIZE LOOP COUNTER */
      if (pw > 0.0)
      {
         cycles = (mix - trim) / (64.66*pw+32*window);
         cycles = 2.0*(double) (int) (cycles/2.0);
      }
      else
      {
//.........这里部分代码省略.........

void pulsesequence()
{
double	Dtau,Ddelta,dosytimecubed,dosyfrq,delcor,
        del    = getval("del"),
        gstab  = getval("gstab"),
	gt1    = getval("gt1"),
	gzlvl1 = getval("gzlvl1"),
	gzlvlhs   = getval("gzlvlhs"),
	hsgt     = getval("hsgt"),
        satpwr = getval("satpwr"),
        satdly = getval("satdly"),
        satfrq = getval("satfrq");
char	delflag[MAXSTR],sspul[MAXSTR],
        alt_grd[MAXSTR],satmode[MAXSTR],lkgate_flg[MAXSTR];

   getstr("delflag",delflag);
   getstr("satmode",satmode);
   getstr("alt_grd",alt_grd);
   getstr("lkgate_flg",lkgate_flg);
   getstr("sspul",sspul);

       //synchronize gradients to srate for probetype='nano'
       //   Preserve gradient "area"
          gt1 = syncGradTime("gt1","gzlvl1",0.5);
          gzlvl1 = syncGradLvl("gt1","gzlvl1",0.5);

   /* In pulse sequence, minimum del=4.0*pw+3*rof1+gt1+2.0*gstab	*/
   if (del < (4*pw+3.0*rof1+gt1+2.0*gstab))
   {  abort_message("Dbppste error: 'del' is less than %g, too short!",
		(4.0*pw+3*rof1+gt1+2.0*gstab));
   }

   Ddelta=gt1;
   Dtau=2.0*pw+rof1+gstab+gt1/2.0;
   dosyfrq = sfrq;
   dosytimecubed=Ddelta*Ddelta*(del-(Ddelta/3.0)-(Dtau/2.0));
   putCmd("makedosyparams(%e,%e)\n",dosytimecubed,dosyfrq);

   /* Check for the validity of the number of transients selected ! 
   if ( (nt != 4) && ((int)nt%16 != 0) )
   {  abort_message("Dbppste error: nt must be 4, 16, 64 or n*16 (n: integer)!");
   }  */

   /* phase cycling calculation */
   mod2(ct,v1); dbl(v1,v1);    hlv(ct,v2);
   mod2(v2,v3); dbl(v3,v3);    hlv(v2,v2);
   mod2(v2,v4); add(v1,v4,v1);			/*   v1    */
   hlv(v2,v2);  add(v2,v3,v4);			/*   v4    */
   hlv(v2,v2);  mod2(v2,v3);   dbl(v3,v5);
   hlv(v2,v2);  mod2(v2,v3);   dbl(v3,v3);	/*   v3    */
   hlv(v2,v2);  mod2(v2,v6);   add(v5,v6,v5);	/*   v5    */
   hlv(v2,v2);  mod2(v2,v2);   dbl(v2,v2);	/*   v2    */

   assign(v1,oph);  dbl(v2,v6);      sub(oph,v6,oph);
   add(v3,oph,oph); sub(oph,v4,oph); dbl(v5,v6);
   add(v6,oph,oph); mod4(oph,oph);		/*receiver phase*/
   add(two,oph,oph);

   mod2(ct,v10);	/* gradients change sign at odd transients */

   if (ni > 1.0)
   {  abort_message("Dbppste is a 2D, not a 3D dosy sequence:  please set ni to 0 or 1");
   }

   /* equilibrium period */
   status(A);

     if (sspul[0]=='y')
       {
         zgradpulse(gzlvlhs,hsgt);
         rgpulse(pw,zero,rof1,rof1);
         zgradpulse(gzlvlhs,hsgt);
       }
      if (satmode[0] == 'y')
       {
        if (d1 - satdly > 0) delay(d1 - satdly);
        obspower(satpwr);
        if (satfrq != tof) obsoffset(satfrq);
        rgpulse(satdly,zero,rof1,rof1);
        if (satfrq != tof) obsoffset(tof);
        obspower(tpwr);
        delay(1.0e-5);
      }
     else
     {  delay(d1); }

if (getflag("wet"))
     wet4(zero,one);

   status(B);
   /* first part of bppdel sequence */
   if (delflag[0]=='y')
   {  if (gt1>0 && gzlvl1>0)
      {  rgpulse(pw, v1, rof1, 0.0);		/* first 90, v1 */
         if (lkgate_flg[0] == 'y') lk_hold();   /* turn lock sampling off */
         if (alt_grd[0] == 'y')                
         {  ifzero(v10);
              zgradpulse(gzlvl1,gt1/2.0);
            elsenz(v10);
              zgradpulse(-1.0*gzlvl1,gt1/2.0);
//.........这里部分代码省略.........

void pulsesequence()
{
   double	base,
                corr,
                presat,
                qlvl;
   char         sspul[MAXSTR];


/* LOAD VARIABLES AND CHECK CONDITIONS */
   presat = getval("presat");
   qlvl = getval("qlvl");
   getstr("sspul", sspul);

   base = 180.0 / qlvl;
   initval(2.0 * qlvl, v5);

   if ((rof1 < 9.9e-6) && (ix == 1))
      fprintf(stdout,"Warning:  ROF1 is less than 10 us\n");

/* STEADY-STATE PHASECYCLING */
/* This section determines if the phase calculations trigger off of (SS - SSCTR)
   or off of CT */

   ifzero(ssctr);
      modn(ct, v5, v10);
      divn(ct, v5, v12);
      mod2(ct, v9);
   elsenz(ssctr);
      sub(ssval, ssctr, v14);	/* v14 = 0,...,ss-1 */
      modn(v14, v5, v10);
      divn(v14, v5, v12);
      mod2(v14, v9);
   endif(ssctr);


/* CALCULATE PHASECYCLE */
/* The phasecycle first performs a (2*Q)-step cycle on the third pulse in order
   to select for MQC.  The phasecycle is then adjusted so that the receiver
   goes +- in an alternating fashion.  Second, the 2-step QIS cycle is added
   in.  Third, a 2-step cycle for axial peak suppression is performed on the
   first pulse. */

   assign(v12, v1);
   mod2(v12, v12);		/* v12=quad. image suppression */
   hlv(v1, v1);
   mod2(v1, v1);
   dbl(v1, v1);
   add(v1, v12, v4);
   add(v12, v1, v1);
   assign(v12, v2);
   assign(v12, v3);
   dbl(v9, v9);
   add(v9, v4, v4);
   assign(v4, oph);
   if (phase1 == 2)
      incr(v1);
   if (phase1 == 3)  /* TPPI */
      add(id2, v1, v1);
/* FAD added for phase=1 or phase=2 */
   if ((phase1 == 1) || (phase1 == 2))
   {
      initval(2.0*(double)(d2_index%2),v13);
      add(v1,v13,v1); add(oph,v13,oph);
   }


/* BEGIN ACTUAL PULSE SEQUENCE CODE */
   if (newtrans)
      obsstepsize(base);

   status(A);
   if (sspul[0] == 'y')
   {
      hsdelay(hst + 0.001);
      rgpulse(pw, v1, 1.0e-6, 1.0e-6);
      hsdelay(hst + 0.001);
   }
   if ((d1 - presat) <= hst)
   {
      rcvroff();
      decon();
      hsdelay(presat);
      decoff();
      delay(1.0e-6);
      rcvron();
   }
   else
   {
      hsdelay(d1 - presat);
      decon();
      rcvroff();
      delay(presat);
      decoff();
      delay(1.0e-6);
      rcvron();
   }
   status(B);
      if (newtrans)
         xmtrphase(v10);      /* hardware digital phaseshift */
//.........这里部分代码省略.........

pulsesequence()
{

  double j1min = getval("j1min"),
         j1max = getval("j1max"),
	 pwxlvl180 = getval("pwxlvl180"),
	 pwx180 = getval("pwx180"),
         gzlvl0 = getval("gzlvl0"),
         gt0 = getval("gt0"),
         gzlvlE = getval("gzlvlE"),
         gtE = getval("gtE"),
         EDratio = getval("EDratio"),
         gstab = getval("gstab"),
	 grad1,
	 grad2,
         tau,
	 tauA,
         tauB,
         taumb;

  char   pwx180ad[MAXSTR]; 
  int	 icosel,
         prgcycle = (int)(getval("prgcycle")+0.5),
	 phase1 = (int)(getval("phase")+0.5);

//synchronize gradients to srate for probetype='nano'
//   Preserve gradient "area"
        gtE = syncGradTime("gtE","gzlvlE",1.0);
        gzlvlE = syncGradLvl("gtE","gzlvlE",1.0);

  getstr("pwx180ad", pwx180ad);
  grad1 = gzlvlE;
  grad2 = -1.0*gzlvlE*(EDratio-1)/(EDratio+1);
  tauA = 1/(2*(j1min + 0.146*(j1max - j1min)));
  tauB = 1/(2*(j1max - 0.146*(j1max - j1min)));
  taumb = 1 / (2 * (getval("jnxh")));
  tau = 1 / (j1min+j1max);

  icosel = 1;

  assign(ct,v17);
  assign(zero,v18);
  assign(zero,v19);

  if (getflag("prgflg") && (satmode[0] == 'y') && (prgcycle > 1.5))
    {
        hlv(ct,v17);
        mod2(ct,v18); dbl(v18,v18);
        if (prgcycle > 2.5)
           {
                hlv(v17,v17);
                hlv(ct,v19); mod2(v19,v19); dbl(v19,v19);
           }
     }

  settable(t1, 1, ph1);
  settable(t2, 1, ph2);
  settable(t3, 2, ph3);
  settable(t4, 1, ph4);
  settable(t5, 4, ph5);
  settable(t6, 4, ph6);

  getelem(t1, v17,v1);
  getelem(t2, v17,v2);
  getelem(t3, v17, v3);
  getelem(t4, v17, v4);
  getelem(t5, v17, v5);
  getelem(t6, v17, oph);

  add(oph,v18,oph);
  add(oph,v19,oph);

/*
  mod2(id2, v10);
  dbl(v10,v10);
*/
  initval(2.0*(double)(((int)(d2*getval("sw1")+0.5)%2)),v10);

  add(v3, v10, v3);
  add(oph, v10, oph);

  if ((phase1 == 2) || (phase1 == 5))
  {
    icosel = -1;
    grad1 = gzlvlE*(EDratio-1)/(EDratio+1);
    grad2 = -1.0*gzlvlE;
  }

  status(A);

   delay(5.0e-5);
   if (getflag("sspul"))
        steadystate();

   if (satmode[0] == 'y')
     {
        if ((d1-satdly) > 0.02)
                delay(d1-satdly);
        else
                delay(0.02);
//.........这里部分代码省略.........

void pulsesequence()
{
   double	   mixN = getval("mixN"),
		   gzlvlC = getval("gzlvlC"),
		   gtC = getval("gtC"),
		   gstab = getval("gstab"),
                   selpwrA = getval("selpwrA"),
                   selpwA = getval("selpwA"),
                   gzlvlA = getval("gzlvlA"),
                   gtA = getval("gtA"),
                   selpwrB = getval("selpwrB"),
                   selpwB = getval("selpwB"),
                   gzlvlB = getval("gzlvlB"),
                   gtB = getval("gtB"),
		   selfrq = getval("selfrq"),
		   zqfpw3 = getval("zqfpw3"),
		   zqfpwr3 = getval("zqfpwr3"),
		   gzlvlzq3 = getval("gzlvlzq3"),
		   mixNcorr,
		   sweeppw = getval("sweeppw"),
		   sweeppwr = getval("sweeppwr");
   char		   sweepshp[MAXSTR],
                   selshapeA[MAXSTR],
                   selshapeB[MAXSTR],
		   zqfpat3[MAXSTR];
//synchronize gradients to srate for probetype='nano'
//   Preserve gradient "area"
        gtC = syncGradTime("gtC","gzlvlC",1.0);
        gzlvlC = syncGradLvl("gtC","gzlvlC",1.0);
        gtA = syncGradTime("gtA","gzlvlA",1.0);
        gzlvlA = syncGradLvl("gtA","gzlvlA",1.0);
        gtB = syncGradTime("gtB","gzlvlB",1.0);
        gzlvlB = syncGradLvl("gtB","gzlvlB",1.0);

   getstr("sweepshp",sweepshp);
   getstr("selshapeA",selshapeA);
   getstr("selshapeB",selshapeB);
   getstr("zqfpat3",zqfpat3);

   mixNcorr=0.0;
   if (getflag("Gzqfilt"))
	mixNcorr=getval("zqfpw3");

   hlv(ct,v1); hlv(v1,v1); hlv(v1,v1); hlv(v1,v1); mod4(v1,v1);
   mod4(ct,v2); add(v1,v2,v2);
   hlv(ct,v3); hlv(v3,v3); mod4(v3,v3); add(v1,v3,v3); dbl(v3,v4);
   dbl(v2,oph); add(oph,v4,oph); add(oph,v1,oph);


/* BEGIN THE ACTUAL PULSE SEQUENCE */
   status(A);
   obspower(tpwr);
   delay(5.0e-5);

   delay(5.0e-5);
   if (getflag("sspul"))
        steadystate();

   if (satmode[0] == 'y')
     {
        if ((d1-satdly) > 0.02)
                delay(d1-satdly);
        else
                delay(0.02);
        satpulse(satdly,zero,rof1,rof1);
     }
   else
        delay(d1);

   if (getflag("wet"))
     wet4(zero,one);

   status(B);
      if (selfrq != tof)
       obsoffset(selfrq);
      rgpulse(pw, v1, rof1, rof1);

      if (selfrq != tof)
        obsoffset(selfrq);

        zgradpulse(gzlvlA,gtA);
        delay(gstab);
        obspower(selpwrA);
        shaped_pulse(selshapeA,selpwA,v2,rof1,rof1);
        obspower(tpwr);
        zgradpulse(gzlvlA,gtA);
        delay(gstab);

      if (selfrq != tof)
        delay(2*OFFSET_DELAY);

        zgradpulse(gzlvlB,gtB);
        delay(gstab);
        obspower(selpwrB);
        shaped_pulse(selshapeB,selpwB,v3,rof1,rof1);
        obspower(tpwr);
        zgradpulse(gzlvlB,gtB);
        delay(gstab);

      if (selfrq != tof)
//.........这里部分代码省略.........