# [neutron-mc] Error in component Source_simple.comp ?

Frederik Zilly zilly at hmi.de
Wed Feb 15 16:01:04 CET 2006

```Hello Peter !

Thank you very much for your quick answer.
Indeed, I must have been blind not seeing the square root.
With this square root everything is alright as your jpg shows.

Sorry for this "false alarm".

Cheers,

Frederik.

Peter Willendrup wrote:

> Hello Frederik,
>
> On Tue, 2006-02-14 at 16:12 +0100, Frederik Zilly wrote:
>
>>Points near the center have greater probility than points with a higher
>
>
> Actually, this is not the case. The square root actually gives zero
> probability for points exactly at the origin; more events must be
> distributed over a circular area element at the circle perimeter  than
> at the center to achieve a uniform distribution.
>
> I wrote up a few lines of matlab code to show this in graphics (you need
> the attached hist2d code to reproduce the attached graphics):
>
> chi=2*pi*rand(1e6,1);
> r=sqrt(rand(1e6,1)); % Source radius assumed one metre
> x=r.*cos(chi);
> y=r.*sin(chi);
> mXY=[x y];
> vXEdge = linspace(-1,1,128);
> vYEdge = linspace(-1,1,128);
> img=hist2d(mXY,vXEdge,vYEdge);
> imagesc(img);
>
> As you can see in the jpg graphic, the result of the above code (which
> is equvalent to the c-code in Source_simple) is completely uniform
> across the source area.
>
> Hope this helps?
>
>
> Cheers,
>
> Peter
>
>
> ------------------------------------------------------------------------
>
> %function mHist = hist2d ([vY, vX], vYEdge, vXEdge)
> %2 Dimensional Histogram
> %Counts number of points in the bins defined by vYEdge, vXEdge.
> %size(vX) == size(vY) == [n,1]
> %size(mHist) == [length(vYEdge) -1, length(vXEdge) -1]
> %
> %EXAMPLE
> %   mYX = rand(100,2);
> %   vXEdge = linspace(0,1,10);
> %   vYEdge = linspace(0,1,20);
> %   mHist2d = hist2d(mYX,vYEdge,vXEdge);
> %
> %   nXBins = length(vXEdge);
> %   nYBins = length(vYEdge);
> %   vXLabel = 0.5*(vXEdge(1:(nXBins-1))+vXEdge(2:nXBins));
> %   vYLabel = 0.5*(vYEdge(1:(nYBins-1))+vYEdge(2:nYBins));
> %   pcolor(vXLabel, vYLabel,mHist2d); colorbar
> function mHist = hist2d (mX, vYEdge, vXEdge)
> [nRow, nCol] = size(mX);
> if nCol < 2
>     error ('mX has less than two columns')
> end
>
> nRow = length (vYEdge)-1;
> nCol = length (vXEdge)-1;
>
> vRow = mX(:,1);
> vCol = mX(:,2);
>
> mHist = zeros(nRow,nCol);
>
> for iRow = 1:nRow
>     rRowLB = vYEdge(iRow);
>     rRowUB = vYEdge(iRow+1);
>
>     [mIdxRow] = find (vRow > rRowLB & vRow <= rRowUB);
>     vColFound = vCol(mIdxRow);
>
>     if (~isempty(vColFound))
>
>
>         vFound = histc (vColFound, vXEdge);
>
>         nFound = (length(vFound)-1);
>
>         if (nFound ~= nCol)
>             [nFound nCol]
>             error ('hist2d error: Size Error')
>         end
>
>         [nRowFound, nColFound] = size (vFound);
>
>         nRowFound = nRowFound - 1;
>         nColFound = nColFound - 1;
>
>         if nRowFound == nCol
>             mHist(iRow, :)= vFound(1:nFound)';
>         elseif nColFound == nCol
>             mHist(iRow, :)= vFound(1:nFound);
>         else
>             error ('hist2d error: Size Error')
>         end
>     end
>
> end
>
>
>
>
> ------------------------------------------------------------------------
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> neutron-mc mailing list
> neutron-mc at risoe.dk
> http://mailman.risoe.dk/mailman/listinfo/neutron-mc

```