src/BWRRelax/relaxRF.cc
void Rrf_2s(super_op& LOp, int rank, gen_op* T1s, gen_op* T2s, matrix* J12)

	// Input		LOp   : Current relaxation superoperator
	//			w     : Vector of system energies
	//				in rad/sec in the LAB frame
	//			rank  : Rank of the two interactions
	// 			T1s   : Spin tensor components, 1st interaction
	// 			T2s   : Spin tensor components, 2nd interaction
	//			J12   : Matrix of reduced spectral density 
	//				function values for 2 interactions
	//				times the two interaction constants
	// Return		void  : Level 2 relax superop for interactions
	//			        1 & 2 added to input superoperator LOp
	//
	//				       LOp    = LOp   + R2
	//				          out	   in     12
	//
	// Note				T1s, T2s, & LOp assumed in proper bases

/*                   --- ---
                     \   \       m    m       -m
             R2    = /   /   (-1)  [ T   , [ T   ,  ] ] J  (w - m*W  )
  	       1,2   --- ---          1,p     2,p        12  p     rf
  		      p   m                                                  */

// where p sums over transitions, m sums over angular momentum components,
// and J values include xi1 and xi2.  The Level 2 treatment assumes that there
// are no degenerate transitions in the system.

// sosi UNDER CONSTRUCTION: THIS  DOESN'T WORK YET
   {
   int hs = T1s[0].dim();			// Get Hilbert space size
   matrix mx1(hs, hs, 0.0);			// Two Zeroed Hilbert space mxs
   matrix mx2(hs, hs, 0.0);
   basis bs = T1s[0].get_basis();		// Current working basis
   gen_op nullOp;				// A NULL operator
//   double Jwp;				// Spectral density, transition p
   complex Ta_aa;
   gen_op *T1ps, *T2ps;				// Components of the two spin tensors
   T1ps = new gen_op[2*rank+1]; 		// (2*l+1 of these for each p)
   T2ps = new gen_op[2*rank+1];			// (2*l+1 of these for each p)
   int m=0, k=0;				// Angular momentum indicies
   for(int a=0; a<hs; a++)			// Sum over transitions a-aa: p
     for(int aa=0; aa<hs; aa++)
       {
       for(m=-rank, k=0; m<=rank; m++, k++)	// Perform pre-summation over m
         {					// & set up spin tensor operator
         T1ps[k] = nullOp;			// components this transition
         Ta_aa = T1s[k].get(a,aa);		// 
         if(Re(Ta_aa) || Im(Ta_aa))		//                       m
           { 					//           T1ps[k] == T
           mx1.put(Ta_aa, a, aa); 		//                       1,p
           T1ps[k] = gen_op(mx1,bs);
           mx1.put(0.0, a, aa);
           }
         T2ps[k] = nullOp;
         Ta_aa = T2s[k].get(a,aa);		//                       m
         if(Re(Ta_aa) || Im(Ta_aa))		//           T2ps[k] == T
           {					//                       2,p
           mx2.put(Ta_aa, aa, a);		//
           T2ps[k] = gen_op(mx2,bs);		// These are sparse or null
           mx2.put(0.0, aa, a);			// so insure null components are
           }					// recognized (& not used).
         }
// Need a vector of Jwp[m] here instead
//Jwp = Re(J12.get(a,aa));			// Get spectral density for this transition
//R_CC_0_trans(T1ps, T2ps, LOp, rank, Jwp);// Now, perform sum over m, and add to LOp
// sosi - the above is still goofy 
       }
   return;
   if(LOp.dim() && J12==NULL) hs=0;		// Compiler likes these used 
   }
// ------------------ Level 2 via Double Commutators --------------------
//                      (No Degenerate Transitions)
Level 2 via Double Commutators
(No Degenerate Transitions)
Level 2 via Double Commutators ( No Degenerate Transitions )
void Rrf_2s(super_op& LOp, int rank, gen_op* T1s, gen_op* T2s, matrix* J12)
					       LOp    = LOp   + R2
					          out	   in     12
	 Note				T1s, T2s, & LOp assumed in proper bases
 where p sums over transitions, m sums over angular momentum components,
 and J values include xi1 and xi2.  The Level 2 treatment assumes that there
 are no degenerate transitions in the system.
 sosi UNDER CONSTRUCTION: THIS  DOESN'T WORK YET
   {
   int hs = T1s[0].dim();			 Get Hilbert space size
   matrix mx1(hs, hs, 0.0);			 Two Zeroed Hilbert space mxs
   matrix mx2(hs, hs, 0.0);
   basis bs = T1s[0].get_basis();		 Current working basis
   gen_op nullOp;				 A NULL operator
   double Jwp;				 Spectral density, transition p
   complex Ta_aa;
   gen_op *T1ps, *T2ps;				 Components of the two spin tensors
   T1ps = new gen_op[2*rank+1]; 		 (2*l+1 of these for each p)
   T2ps = new gen_op[2*rank+1];			 (2*l+1 of these for each p)
   int m=0, k=0;				 Angular momentum indicies
   for(int a=0; a<hs; a++)			 Sum over transitions a-aa: p
     for(int aa=0; aa<hs; aa++)
       {
       for(m=-rank, k=0; m<=rank; m++, k++)	 Perform pre-summation over m
         {					 & set up spin tensor operator
         T1ps[k] = nullOp;			 components this transition
         Ta_aa = T1s[k].get(a,aa);		 
         if(Re(Ta_aa) || Im(Ta_aa))		                       m
           { 					           T1ps[k] == T
           mx1.put(Ta_aa, a, aa); 		                       1,p
           T1ps[k] = gen_op(mx1,bs);
           mx1.put(0.0, a, aa);
           }
         T2ps[k] = nullOp;
         Ta_aa = T2s[k].get(a,aa);		                       m
         if(Re(Ta_aa) || Im(Ta_aa))		           T2ps[k] == T
           {					                       2,p
           mx2.put(Ta_aa, aa, a);		
           T2ps[k] = gen_op(mx2,bs);		 These are sparse or null
           mx2.put(0.0, aa, a);			 so insure null components are
           }					 recognized (& not used).
         }
 Need a vector of Jwp[m] here instead
Jwp = Re(J12.get(a,aa));			 Get spectral density for this transition
R_CC_0_trans(T1ps, T2ps, LOp, rank, Jwp); Now, perform sum over m, and add to LOp
 sosi - the above is still goofy 
       }
   return;
   if(LOp.dim() && J12==NULL) hs=0;		 Compiler likes these used 
   }
void Rrf_2s ( super_op & LOp , int rank , gen_op * T1s , gen_op * T2s , matrix * J12 ) { int hs = T1s [ 0 ] . dim ( ) ; matrix mx1 ( hs , hs , 0.0 ) ; matrix mx2 ( hs , hs , 0.0 ) ; basis bs = T1s [ 0 ] . get_basis ( ) ; gen_op nullOp ; complex Ta_aa ; gen_op * T1ps , * T2ps ; T1ps = new gen_op [ 2 * rank + 1 ] ; T2ps = new gen_op [ 2 * rank + 1 ] ; int m = 0 , k = 0 ; for ( int a = 0 ; a < hs ; a ++ ) for ( int aa = 0 ; aa < hs ; aa ++ ) { for ( m = - rank , k = 0 ; m <= rank ; m ++ , k ++ ) { T1ps [ k ] = nullOp ; Ta_aa = T1s [ k ] . get ( a , aa ) ; if ( Re ( Ta_aa ) || Im ( Ta_aa ) ) { mx1 . put ( Ta_aa , a , aa ) ; T1ps [ k ] = gen_op ( mx1 , bs ) ; mx1 . put ( 0.0 , a , aa ) ; } T2ps [ k ] = nullOp ; Ta_aa = T2s [ k ] . get ( a , aa ) ; if ( Re ( Ta_aa ) || Im ( Ta_aa ) ) { mx2 . put ( Ta_aa , aa , a ) ; T2ps [ k ] = gen_op ( mx2 , bs ) ; mx2 . put ( 0.0 , aa , a ) ; } } } return ; if ( LOp . dim ( ) && J12 == NULL ) hs = 0 ; }
Level 2 via Double Commutators
(No Degenerate Transitions)
Level 2 via Double Commutators ( No Degenerate Transitions )
// Input		LOp   : Current relaxation superoperator //			w     : Vector of system energies //				in rad/sec in the LAB frame //			rank  : Rank of the two interactions // 			T1s   : Spin tensor components, 1st interaction // 			T2s   : Spin tensor components, 2nd interaction //			J12   : Matrix of reduced spectral density  //				function values for 2 interactions //				times the two interaction constants // Return		void  : Level 2 relax superop for interactions //			        1 & 2 added to input superoperator LOp // //				       LOp    = LOp   + R2 //				          out	   in     12 // // Note				T1s, T2s, & LOp assumed in proper bases /*                   --- ---
                     \   \       m    m       -m
             R2    = /   /   (-1)  [ T   , [ T   ,  ] ] J  (w - m*W  )
  	       1,2   --- ---          1,p     2,p        12  p     rf
  		      p   m                                                  */ // where p sums over transitions, m sums over angular momentum components, // and J values include xi1 and xi2.  The Level 2 treatment assumes that there // are no degenerate transitions in the system. // sosi UNDER CONSTRUCTION: THIS  DOESN'T WORK YET // Get Hilbert space size // Two Zeroed Hilbert space mxs // Current working basis // A NULL operator //   double Jwp;				// Spectral density, transition p // Components of the two spin tensors // (2*l+1 of these for each p) // (2*l+1 of these for each p) // Angular momentum indicies // Sum over transitions a-aa: p // Perform pre-summation over m // & set up spin tensor operator // components this transition //  //                       m //           T1ps[k] == T //                       1,p //                       m //           T2ps[k] == T //                       2,p // // These are sparse or null // so insure null components are // recognized (& not used). // Need a vector of Jwp[m] here instead //Jwp = Re(J12.get(a,aa));			// Get spectral density for this transition //R_CC_0_trans(T1ps, T2ps, LOp, rank, Jwp);// Now, perform sum over m, and add to LOp // sosi - the above is still goofy  // Compiler likes these used 