anchor.diff:             double precision, dimension(size(ncgs),11) :: ktemp
              double precision, dimension(size(ncgs),11),intent(out) :: ktotal
   !           ones = 1d0
  +!           write(6,*) 'maxjutt: ',size(maxjutt_args),size(weights_arr)
              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
  +!           write(6,*) 'maxjutt weights: ',weights_arr,delta
              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)
  
positive.code: !    ktemp  = polsynchth(nu,weights_arr(i)*ncgs,bcgs,thetamin*delta_arr(i),e%incang,Kth)           
  !    ktotal = ktotal + ktemp
  !return K to emis.f90
  
             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
           end subroutine calc_maxjutt_subroutine
  
           end module calc_maxjutt
  
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
  