Sorry again: This ofc is wrong: volflow(eta,0.05)=2*broadcast(+, eta.^2-eta.^4/2 , sum( broadcast(*, broadcast(/,2*eta,r.^2).*broadcast(-,broadcast(*,2./r,broadcast(/,besselj1(broadcast(*,r,eta)),besselj0(r))),eta) , exp(-broadcast(*,r.^2,x_s))) , 3) )
correct: volflow(eta)=vol_flow(eta,0.05)
