function varrtp=rtpvariable(magc,incv);
% *****************************************************************
% *** RTP with variable geomagnetic parameters
% *** Uses a Taylor series expansion in the space domain
% *** GRJ Cooper July 2002
% ******************************************************************
% *** magc  ; input data (even no of datapoints in each dir)
% *** incv  ; dataset of field inc over dataset - must be exactly
% *** the same size as magc. If omitted then inc is specified only
% *** at the dataset corners, and interpolated from there
% ******************************************************************
% *** Type rtpvariable with no input parameters for demo application
% ******************************************************************

dtr=pi/180; azimuth=0;

% interpolate inc variation over dataset using linear polynomial

if nargin==0
 disp('Loading demo data');
 load rtpdemodata;
end;
[nr,nc]=size(magc); 
if nargin<2 % Specify inc variation only at corners
 disp('Generating grid of inclination data from values at grid corners');
 a=zeros(4,3); yval=zeros(4,1);
 a(1,1)=1; a(1,2)=nr; a(1,3)=1; yval(1,1)=-30;     
 a(2,1)=1; a(2,2)=nr; a(2,3)=nc; yval(2,1)=-30;  
 a(3,1)=1; a(3,2)=1; a(3,3)=1; yval(3,1)=-90;    
 a(4,1)=1; a(4,2)=1; a(4,3)=nc; yval(4,1)=-90; 
 p=pinv(a)*yval;
 [xi,yi]=meshgrid(1:nc,1:nr); 
 incv=p(1)+p(2)*yi+p(3)*xi;
end;
incv=incv*dtr; inc=incv(nr/2,nc/2); % inclination at centre of grid

nmax=max([nr nc]); npts=(2^nextpow2(nmax));
cdiff=floor((npts-nc)/2); rdiff=floor((npts-nr)/2);
disp('Computing rtp');
rtpdata=rtp(magc,npts,nr,nc,rdiff,cdiff,azimuth,inc);
incv=incv-inc;
disp('Computing dx1');
drtp1=taylortp1(magc,npts,nr,nc,rdiff,cdiff,azimuth,inc);
disp('Computing dx2');
drtp2=taylortp1(drtp1,npts,nr,nc,rdiff,cdiff,azimuth,inc);
disp('Computing dx3');
drtp3=taylortp1(drtp2,npts,nr,nc,rdiff,cdiff,azimuth,inc);
disp('Computing rtpvar');
varrtp=rtpdata+drtp1.*incv+0.5*drtp2.*incv.^2+0.1666666*drtp3.*incv.^3; % 3 terms - modify as required
figure(1); 
subplot(1,2,1); imagesc(magc); axis equal; axis tight; axis xy; title('Data'); axis off; 
hold on; contour(magc,6,'k'); hold off;
subplot(1,2,2); imagesc((incv+inc)/dtr); axis equal; axis tight; axis xy; title('Inclination Variation');
axis off; colorbar;
currfig=get(0,'CurrentFigure'); set(currfig,'numbertitle','off');
set(currfig,'name','RTP Datasets'); 
figure(2);
subplot(1,3,1); imagesc(drtp1); axis equal; axis tight; axis xy; title('1st Derivative'); axis off; 
subplot(1,3,2); imagesc(drtp2); axis equal; axis tight; axis xy; title('2nd Derivative'); axis off; 
subplot(1,3,3); imagesc(drtp3); axis equal; axis tight; axis xy; title('3rd Derivative'); axis off; 
currfig=get(0,'CurrentFigure'); set(currfig,'numbertitle','off');
set(currfig,'name','RTP Derivatives'); 
figure(3);
subplot(1,2,1); imagesc(rtpdata); axis equal; axis tight; axis xy; title('After Standard RTP'); axis off; 
hold on; contour(rtpdata,6,'k'); hold off;
subplot(1,2,2); imagesc(varrtp); axis equal; axis tight; axis xy; title('c) After Differential RTP'); axis off;
hold on; contour(varrtp,6,'k'); hold off;
colormap('jet'); 
currfig=get(0,'CurrentFigure'); set(currfig,'numbertitle','off');
set(currfig,'name','RTP with Variable Inclination Over Survey Area'); 

%**************************************************************************
function drtp=taylortp1(magc,npts,nr,nc,rdiff,cdiff,azimuth,inc); % Rtp derivatives

dtheta=0.01; dtr=pi/180;
rtpdata1=rtp(magc,npts,nr,nc,rdiff,cdiff,azimuth,inc);
rtpdata2=rtp(magc,npts,nr,nc,rdiff,cdiff,azimuth,inc+dtheta);
drtp=(rtpdata2-rtpdata1)/dtheta;

% *************************************************************************
function rtpdata=rtp(magc,npts,nr,nc,rdiff,cdiff,azimuth,inc);

gf=taper2d(magc,npts,nc,nr,cdiff,rdiff);
f=fft2(gf,npts,npts);
fx=f;

% Direction cosines

L=cos(azimuth)*cos(inc);
M=sin(azimuth)*cos(inc);
N=sin(inc);

wn=2.0*pi/(npts-1); 
for I=1:npts/2
 for J=1:npts/2
  U=I*wn; V=J*wn;
  temp=U*U+V*V;
  temp1=L*U+M*V;
  temp2=N*N*temp-temp1*temp1;
  base=temp2*temp2+4.0*N*N*temp1*temp1*temp;
  Rpart=temp*temp2/base;
  Ipart=-2.0*N*temp1*(temp^1.5)/base;
  temp3=real(f(I,J));
  temp4=imag(f(I,J));
  fx(I,J)=real(Rpart*temp3-temp4*Ipart)+i*real(Rpart*temp4+Ipart*temp3);
  
  temp1=-L*U+M*V;
  temp2=N*N*temp-temp1*temp1;
  base=temp2*temp2+4.0*N*N*temp1*temp1*temp;
  Rpart=temp*temp2/base;
  Ipart=-2.0*N*temp1*(temp^1.5)/base;
  temp3=real(f(npts-I+1,J));
  temp4=imag(f(npts-I+1,J));
  fx(npts-I+1,J)=real(Rpart*temp3-temp4*Ipart)+i*real(Rpart*temp4+Ipart*temp3);

  temp1=L*U-M*V;
  temp2=N*N*temp-temp1*temp1;
  base=temp2*temp2+4.0*N*N*temp1*temp1*temp;
  Rpart=temp*temp2/base;
  Ipart=-2.0*N*temp1*(temp^1.5)/base;
  temp3=real(f(I,npts-J+1));
  temp4=imag(f(I,npts-J+1));
  fx(I,npts-J+1)=real(Rpart*temp3-temp4*Ipart)+i*real(Rpart*temp4+Ipart*temp3);

  temp1=-L*U-M*V;
  temp2=N*N*temp-temp1*temp1;
  base=temp2*temp2+4.0*N*N*temp1*temp1*temp;
  Rpart=temp*temp2/base;
  Ipart=-2.0*N*temp1*(temp^1.5)/base;
  temp3=real(f(npts-I+1,npts-J+1));
  temp4=imag(f(npts-I+1,npts-J+1));
  fx(npts-I+1,npts-J+1)=real(Rpart*temp3-temp4*Ipart)+i*real(Rpart*temp4+Ipart*temp3);
 end;
end;
fxinv=ifft2(fx); rtpdata=real(fxinv(1+rdiff:nr+rdiff,1+cdiff:nc+cdiff));

%**************************************************************************
function gf=taper2d(g,npts,nc,nr,cdiff,rdiff); 

gf=zeros(npts); gf(rdiff+1:rdiff+nr,cdiff+1:cdiff+nc)=g;
 for J=cdiff+1:cdiff+nc
  for I=1:rdiff   
   gf(I,J)=gf(2*rdiff-I+1,J)*((1+sin(-pi/2+(I-1)*pi/rdiff))*0.5); 
   gf(npts-I+1,J)=gf(npts-2*rdiff+I,J)*((1+sin(-pi/2+(I-1)*pi/rdiff))*0.5);                           
  end;
 end;
 
 for I=rdiff+1:rdiff+nr
  for J=1:cdiff   
   gf(I,J)=gf(I,2*cdiff-J+1)*((1+sin(-pi/2+(J-1)*pi/(cdiff)))*0.5);     
   gf(I,npts-J)=gf(I,npts-2*cdiff+J)*((1+sin(-pi/2+(J-1)*pi/(cdiff)))*0.5);                                                     
  end;
 end;

 for J=1:cdiff
  for I=1:rdiff
   gf(I,J)=gf(rdiff+1,J)*((1+sin(-pi/2+(I-1)*pi/(rdiff)))*0.5); 
   gf(npts-I,J)=gf(npts-rdiff,J)*((1+sin(-pi/2+(I-1)*pi/(rdiff)))*0.5); 
  end;
 end;

 for J=cdiff+nc:npts
  for I=1:rdiff
   gf(I,J)=gf(rdiff+1,J)*((1+sin(-pi/2+(I-1)*pi/(rdiff)))*0.5); 
   gf(npts-I,J)=gf(npts-rdiff,J)*((1+sin(-pi/2+(I-1)*pi/(rdiff)))*0.5); 
  end;
 end;