anchor.diff:  	write (*,*) 'case 1'
   	do p_counter = 1,size(propagators,4)
   		do t_counter = 1,size(propagators,3)
  -			call sandwich_internal(&
  -				rho(:,:,t_counter,p_counter),&
  -				propagators(:,:,t_counter,p_counter),&
  -				propagators(:,:,t_counter,p_counter))
  +			call sandwich_internal(rho(:,:,t_counter,p_counter),propagators(:,:,t_counter,p_counter),propagators(:,:,t_counter,p_counter))
   		enddo
   	enddo
   elseif (size(rho,4).eq.1.and.size(rho,3).eq.1) then
  
positive.code: 	do p_counter = 1,size(propagators,4)
  		do t_counter = 1,size(propagators,3)
  			!write (*,*) 'multiply',t_counter,p_counter
  			propagators(:,:,t_counter,p_counter) = matmul(&
  				propagators(:,:,t_counter,p_counter),&
  				matmul(&
  					rho(:,:,1,1) , conjg(transpose(&
  						propagators(:,:,t_counter,p_counter)&
  						))&
  					))
  		enddo
  	enddo
  else
  
negative.code: 		do t_counter = 1,size(time_dependence,2)
  			temp_exp = 0
  			do op_counter = 1,size(time_dependence,1)
  				temp_exp = temp_exp + p(op_counter,p_counter) * time_dependence(op_counter,t_counter) * operators(:,:,op_counter)
  			end do
  			call complexexp(temp_exp) ! this passes a HERMITIAN argument, and gets back the exponent of the 2*pi*i* that argument
  			! for now, i'm not actually propagating, just exponentiating
  			if (t_counter.gt.1) then
  				propagators(:,:,1,p_counter) = matmul( temp_exp, propagators(:,:,1,p_counter))
  			else
  				propagators(:,:,1,p_counter) = temp_exp
  			endif
  