src/Level2/TrnsTable1D.cc
/* This routine phase corrects all the transitions in matrix mx using a
   0th order and 1st order phase correction.  The 0th order finds the
   transition closest to w0, takes the phase of that transition, and then
   removes that phase component from all transitions.  The 1st order finds
   the transition closest to w1 and looks at its phase.  It then assumes that
   there is some constant tp such that

                             (w1-w0)*tp = phase(w1) 

   and calculates tp, noting that and addition of 360 degrees to the phase will
   suffice, i.e. that there is another factor, ntp, where

                               (w1-w0)*ntp = 360

   such that

                      (w1-w0)*(tp+n*ntp) = phase(w1) + n*360 

   where n is an integer.  Any value of n will produce a proper phase value
   of (w1-w0)*(tp+n*ntp) that would undo the phase at w1.  The only requirement
   is to select which value of n to use.  To to that, this routine picks the
   biggest transition, say Wtr, and sees which value of n minimizes

                      (Wtr-w0)*(tp+n*ntp) - phase(Wtr) -----> 0

   where n is kept between the range [0, order].  Once found, all transitions
   are phase corrected by adjusting their complex intensities according to

  		    Inew(W) = Iold(W) * exp[-i*(W-w0)*(tp+n*ntp)]            */
This routine phase corrects all the transitions in matrix mx using a
0th order and 1st order phase correction.  The 0th order finds the
transition closest to w0, takes the phase of that transition, and then
removes that phase component from all transitions.  The 1st order finds
the transition closest to w1 and looks at its phase.  It then assumes that
there is some constant tp such that



and calculates tp, noting that and addition of 360 degrees to the phase will
suffice, i.e. that there is another factor, ntp, where



such that










This routine phase corrects all the transitions in matrix mx using a 0th order and 1st order phase correction . The 0th order finds the transition closest to w0 takes the phase of that transition and then removes that phase component from all transitions . The 1st order finds the transition closest to w1 and looks at its phase . It then assumes that there is some constant tp such that and calculates tp noting that and addition of 360 degrees to the phase will suffice i . e . that there is another factor ntp where such that
complex TTable1D::pcorrect(double& w0, double w1, int order)
  {
  double pi2 = 2.0*PI;				
  double w0rad = w0*pi2;			
  double w1rad = w1*pi2;			
  double Wtr;					
  double Wdel0, Wdelmax0 = HUGE_VAL;		
  double Wdel1, Wdelmax1 = HUGE_VAL;		
  int i, i0max=0, i1max=0;
  for(i=0; i<rows(); i++)			
    {
    Wtr = getIm(i,0);				
    Wdel0 = w0rad - Wtr;			
    if(fabs(Wdel0) < Wdelmax0)			
       {					
       Wdelmax0 = fabs(Wdel0);			
       i0max = i;
       }
    Wdel1 = w1rad - Wtr; 			
    if(fabs(Wdel1) < Wdelmax1)			
       {					
       Wdelmax1 = fabs(Wdel1);			
       i1max = i;
       }
    }
  w0rad = getIm(i0max, 0);			
  if(i0max == i1max)				
    {						
    i1max = i0max-1;
    if(i1max < 0) i1max = i0max+1;
    }
  w1rad = getIm(i1max, 0);			
  complex z = get(i0max,1);			
  double phase0 = phase(z);			
  complex pc = norm(z)/z;			 correction
  for(i=0; i<rows(); i++)			
    {
    z = get(i,1);				
    z *= pc;					
    put(z,i,1);					
    }
  z = get(i1max,1);				
  double phase1 = phase(z);			
  if(phase1 < 0) phase1 += pi2;			
  double tp = phase1/(w1rad-w0rad);		
  double ntp = pi2/(w1rad-w0rad);		
  int maxi=0;					
  double Imax = -HUGE_VAL;				
  for(i=0; i<rows(); i++)			
    {
    if(i!=i0max && i!=i1max)			
      {
      if(norm(get(i,1)) > Imax)			
        {					
        Imax = norm(get(i,1));
        maxi = i;
        }
      }
    }
  double wrad = getIm(maxi,0);			
  double prad = phase(get(maxi,1));		
  if(prad < 0) prad += pi2;			
  double delp, pp, pmin=HUGE_VAL;
  int nmin = 0;
  for(i=0; i<order; i++)			
    {
    pp = (wrad-w0rad)*(tp+i*ntp);		
    delp = fabs(prad-pp); 			
    if(delp < pmin)				
      {						
      pmin = delp;				
      nmin = i;
      }
    }
  double tconst = tp+(double(nmin)*ntp);
  double delW;
  for(i=0; i<rows(); i++)			
    {
    Wtr = getIm(i,0);				
    delW = Wtr - w0rad;				
    pp = delW*(tp+nmin*ntp);			
    z = complex(0, pp);				
    z = z.Zexp();				
    pc = norm(z)/z;				
    z = get(i,1);				
    z *= pc;					
    put(z,i,1);					
    }
  w0 = w0rad;					
  return complex(phase0, tconst);		
  }
complex TTable1D :: pcorrect ( double & w0 , double w1 , int order ) { double pi2 = 2.0 * PI ; double w0rad = w0 * pi2 ; double w1rad = w1 * pi2 ; double Wtr ; double Wdel0 , Wdelmax0 = HUGE_VAL ; double Wdel1 , Wdelmax1 = HUGE_VAL ; int i , i0max = 0 , i1max = 0 ; for ( i = 0 ; i < rows ( ) ; i ++ ) { Wtr = getIm ( i , 0 ) ; Wdel0 = w0rad - Wtr ; if ( fabs ( Wdel0 ) < Wdelmax0 ) { Wdelmax0 = fabs ( Wdel0 ) ; i0max = i ; } Wdel1 = w1rad - Wtr ; if ( fabs ( Wdel1 ) < Wdelmax1 ) { Wdelmax1 = fabs ( Wdel1 ) ; i1max = i ; } } w0rad = getIm ( i0max , 0 ) ; if ( i0max == i1max ) { i1max = i0max - 1 ; if ( i1max < 0 ) i1max = i0max + 1 ; } w1rad = getIm ( i1max , 0 ) ; complex z = get ( i0max , 1 ) ; double phase0 = phase ( z ) ; complex pc = norm ( z ) / z ; for ( i = 0 ; i < rows ( ) ; i ++ ) { z = get ( i , 1 ) ; z *= pc ; put ( z , i , 1 ) ; } z = get ( i1max , 1 ) ; double phase1 = phase ( z ) ; if ( phase1 < 0 ) phase1 += pi2 ; double tp = phase1 / ( w1rad - w0rad ) ; double ntp = pi2 / ( w1rad - w0rad ) ; int maxi = 0 ; double Imax = - HUGE_VAL ; for ( i = 0 ; i < rows ( ) ; i ++ ) { if ( i != i0max && i != i1max ) { if ( norm ( get ( i , 1 ) ) > Imax ) { Imax = norm ( get ( i , 1 ) ) ; maxi = i ; } } } double wrad = getIm ( maxi , 0 ) ; double prad = phase ( get ( maxi , 1 ) ) ; if ( prad < 0 ) prad += pi2 ; double delp , pp , pmin = HUGE_VAL ; int nmin = 0 ; for ( i = 0 ; i < order ; i ++ ) { pp = ( wrad - w0rad ) * ( tp + i * ntp ) ; delp = fabs ( prad - pp ) ; if ( delp < pmin ) { pmin = delp ; nmin = i ; } } double tconst = tp + ( double ( nmin ) * ntp ) ; double delW ; for ( i = 0 ; i < rows ( ) ; i ++ ) { Wtr = getIm ( i , 0 ) ; delW = Wtr - w0rad ; pp = delW * ( tp + nmin * ntp ) ; z = complex ( 0 , pp ) ; z = z . Zexp ( ) ; pc = norm ( z ) / z ; z = get ( i , 1 ) ; z *= pc ; put ( z , i , 1 ) ; } w0 = w0rad ; return complex ( phase0 , tconst ) ; }
This routine phase corrects all the transitions in matrix mx using a
0th order and 1st order phase correction.
This routine phase corrects all the transitions in matrix mx using a 0th order and 1st order phase correction .
//		     First Nearest Transitions To w0 and w1 // Required constant // Zero pivot in radians // 1st pivot in radians // Transition frequency // Use to look for w0 // Use to look for w1 // Loop rows // 	Transition i freq. //	Freq. Difference //	Look for a 0th order //	pivot frequency close //	by w0rad	 //	Freq. Difference //	Look for 1st order //	pivot frequency close //	by w1rad // Set zero pivot transition // Insure both pivots aren't // the same point // Set 1st pivot transition //		   Now Apply The Zeroth Order Phase Correction //	         (Transition w0 Will Have Zero Phase Hereafter) // Zero order intensity // Zero order phase // Zero order phase correction // Loop rows // 	Transition intensity // 	Adjust transition phase // 	Reset the intensity //		    Set Up For First Order Phase Correction // 1st order intensity // 1st order phase // Keep phase positive // Base time constant for phase // Adjustable time constant //        Find Largest Transition Which is Not One Of The Two Pivots //	      (We'll Deal With Case Of <= Two Transitions Later) // Index of large transition // Largest transition intensity // Loop transitions // 	Don't use a pivot //	Check intensity //	and store if large // Frequency of large transition // Phase of large transition // Keep phase positive //                 Find Best Value of "n" For 1st Order Phasing // Loop through allowed orders //	Phase made at wrad, n=i //	See if good adjustment //	If its good, store this  //	value of i & use it if //	nothing better is found //		     Apply The First Order Phase Correction // Loop rows // 	Transition i freq. //	Freq. Difference // 	Needed phase correction //	-i*phi //	z = exp(-i*phi) //	Phase correction factor // 	Transition intensity // 	Adjust transition phase // 	Reset (phase) intensity // Set w0 to be pivot point // Return phase factors