Euler rotation confusion

Hello,
I am getting very confused.
I have a list of supernovae with galactic coordinates l (from 0->360) and b (from -180->180).
I am trying to write a program in Matlab that rotates the coordinates by 2 Euler angles. Firstly I change l and b into X,Y,Z i.e change from spherical to cartesian coordinate system.
Then I have the rotation matrix:

[X',Y',Z']=[X,Y,Z]* [cos(a) -sin(a) 0 ; -sin(a) cos(a) 0; 0 0 1]*[1 0 0 ; 0 cos(c) -sin(c) ; 0 sin(c) cos(c)]

I want to return a value (delta) which is a function of the SN properties at each rotation where a ranges from (0->360) in steps of 10 degrees and c ranges from (-90->90) in steps of 5 degrees.

The rotation matrix above only seems to work when the range is set as a from (-180->180) and c from (0->180). This means that I have to change l and b so that b ranges from (0->180) ie b=90-b, and l is in the range (-180->180) so l=l-180.

Then I need to plot the rotation values (that are colour coordinated according to the returned value delta) onto a plot that ranges a from (-180->180) and c that ranges from (-90->90). So I need to change c so that
c=90-c.

However, I think that I am doing something wrong somewhere.
I was wondering if anyone had done something similar to this before or if they can see what the problem is, OR alternatively if anyone knows the Euler rotation matrix that will have the range of a(0->360) and c(-90->90) as I have searched but can't find anything.