function f=calbaz(stlat,stlon,slat,slon)

% function to calculate back azimuth
% stlat; the latitude of the station in degree
% stlon; the longitude of the station in degree
% slat; the latitude of the source in degree
% slat; the longitude of the source in degree

aa=pi/180;
b=90-slat;
c=90-stlat;
if stlon < 0
   stlon=360+stlon;
end   
if slon < 0 
   slon=360+slon;
end   

tbig1=max(stlon,slon);
tbig2=min(stlon,slon);
tempb=tbig1-tbig2;
biga=tempb;
if tempb >+ 180 
   biga=360.-tbig1+tbig2;
end   

	
dd=cos(b*aa)*cos(c*aa)+sin(b*aa)*sin(c*aa)*cos(biga*aa);
a=acos(dd)/aa;
dd1=(cos(b*aa)-cos(a*aa)*cos(c*aa))/(sin(a*aa)*sin(c*aa));
bazi1=acos(dd1)/aa;
bazi2=360-bazi1;	

	
if stlon <= slon 
   if tempb <= 180
   	f=bazi1;
	else
	   f=bazi2;
   end

else
        
	if tempb <= 180
	   f=bazi2;
	else
	   f=bazi1;
   end
end
