function [potini] = cal_potini(tevla,tevlo,refbaz)

%CAL_POTINI: function which is generating potential initial position  
% with respect to reference back azimuth.
% algorithm: 1. Transform given location to north pole.
%            2. Apply pre-defined position from North pole.
%            3. Transform back to original coordinate.
%
% Usage: [potini]=cal_potini(tevla,tevlo,refbaz)

rat1=pi/180;
rat2=180/pi;
baz=10;
epi=10;

% Position of assumed location of the event in Cartesian coordinate 
% under assuming radius of 1. xtev,ytev,ztev 

[xtev,ytev,ztev]=sph2cart(tevla*rat1,tevlo*rat1,1);
xx1=[xtev,ytev,ztev]';

% transforming the position of event to north pole.
azel = [tevlo-90 0];
alpha = 90-tevla;

% Position of event in Cartesian coordinate under assuming station is at NP.

k = 1;
tbaz = refbaz-50:10:refbaz+50;
len_baz = length(tbaz);

for i = 1 : len_baz
for j = 10:10:180
   [xteve,yteve,zteve]=sph2cart((360-tbaz(i))*rat1,(90-j)*rat1,1);
%   [xteve,yteve,zteve]=sph2cart((tevla+tbaz(i))*rat1,(90-j)*rat1,1);
   xxteve(:,k)=[xteve;yteve;zteve];
   k=k+1;
end
end

nmax=k-1;

% find unit vector for axis of rotation
theta = pi*azel(1)/180;
phi = pi*azel(2)/180;
u = [cos(phi)*cos(theta); cos(phi)*sin(theta); sin(phi)];

alph = alpha*pi/180;
cosa = cos(alph);
sina = sin(alph);
vera = 1 - cosa;
x = u(1);
y = u(2);
z = u(3);
xrot = [cosa+x^2*vera x*y*vera-z*sina x*z*vera+y*sina; ...
       x*y*vera+z*sina cosa+y^2*vera y*z*vera-x*sina; ...
       x*z*vera-y*sina y*z*vera+x*sina cosa+z^2*vera]';

xx2=(xrot*xxteve);
[th,phi,r]=cart2sph(xx2(1,:),xx2(2,:),xx2(3,:));

evlo=th*rat2;
evla=phi*rat2;
potini=[evla',evlo'];
