Great Circle Arc Distance Error with Dimensions

I am using this function called grcdistance, which calculates the the great circle arc distance between two points (the event and the station). The function is defined below:

function d=grcdistance(lon1,lat1,lon2,lat2)

% D=GRCAZIM(LON1,LAT1,LON2,LAT2)
% D is the angular distance between (LON1,LAT1) and (LON2,LAT2)
d=acos(sin(lat1).*sin(lat2)+cos(lat1).*cos(lat2).*cos(lon1-lon2));

Now, the problem is that I have ~300 events and 20 stations. When this function tries to calculate the great circle arc distance, I get the following error when I run the script (getparam.m) which calls the grcdistance.m function:

Error using  - 
Matrix dimensions must agree.
Error in grcdistance (line 6)
d=acos(sin(lat1).*sin(lat2)+cos(lat1).*cos(lat2).*cos(lon1-lon2));
Error in getparam (line 15)
dist=grcdistance(evlon,evlat,stlon,stlat)/pi*180;

I know why the error occurs: it's because it's trying to match up the lon1 (of the events, of which there are 300), with the lon2 of the stations (of which there are 20), and there is a mismatch with the matrix. I was wondering what would be the best way to get around this matrix mismatch? For every event (lon1 and lat1), I want to calculate the great circle distance between the each of my 20 stations, so that the output would be a great circle distance for each event-station pair. Thank you for your help.

Answer 1

[James Tursa]

Getting a matrix result using ndgrid:

function d = grcdistance(lon1,lat1,lon2,lat2)
% D = GRCAZIM(LON1,LAT1,LON2,LAT2)
% D is the angular distance between (LON1,LAT1) and (LON2,LAT2)
[LON1 LON2] = ndgrid(lon1,lon2);
[LAT1 LAT2] = ndgrid(lat1,lat2);
d = acos(sin(LAT1).*sin(LAT2)+cos(LAT1).*cos(LAT2).*cos(LON1-LON2));