anchor.diff:  !    ktemp  = polsynchth(nu,weights_arr(i)*ncgs,bcgs,thetamin*delta_arr(i),e%incang,Kth)           
   !    ktotal = ktotal + ktemp
   !return K to emis.f90
  -
  +           ktotal = 0d0
              do i=1,size(weights_arr)
                 ktemp = 0d0
                 call polsynchth(nu,weights_arr(i)*ncgs,bcgs,tcgs_min*delta_arr(i),incang,ktemp)
                 ktotal = ktotal + ktemp
              enddo
  +!           write(6,*) 'maxjutt loop: ',tcgs_min,delta_arr,i,size(weights_arr)
            end subroutine calc_maxjutt_subroutine
   
            end module calc_maxjutt
  
positive.code:            double precision, dimension(size(ncgs),11) :: ktemp
             double precision, dimension(size(ncgs),11),intent(out) :: ktotal
  !           ones = 1d0
             a = 3d0
  !Note derivation assumes relativistic a = 3.0
             delta = maxjutt_args(1) !first argument is delta
             weights_arr = maxjutt_args(2:)/SUM(maxjutt_args(2:)) !normalize
             i_arr = (/ (I, I=0,size(weights_arr)-1) /) !indices
             delta_arr = delta**i_arr !all deltas to be multiplied by thetamin
             dwsum = SUM(weights_arr * delta_arr)
  
negative.code: !6 weights, 1 delta
  !           double precision, dimension(:), intent(in) :: usim
  
             double precision, dimension(size(maxjutt_args) - 1) :: weights_arr, i_arr, delta_arr
             double precision, dimension(size(ncgs)) :: tcgs_min
  !           double precision, dimension(size(ncgs)) :: theta_min
  