%
%  Matlab m-file to load and plot magnetics timescale and magnetics profile
%
%  The file timescale.dat is loaded.
%  The first column of the file is age (in my)
%  The second column of the file is normal (+1) or reversed (-1) field polarity
%
%  The file NBP9707.xydm is loaded.
%  The first column of the file is longitude in degrees.
%  The second column is latitude in degrees.
%  The third column is distance from ridge axis in kilometers.
%  The fourth column is magnetic anomaly in nanotesla.
%
   clear
   clf
   load timescale.dat
   t =timescale(:,1);
   p =timescale(:,2);
%
% plot timescale
%
   plot(t,p)
   axis([0.,10.,-1.5,1.5]);
   xlabel(' Million Years Before Present')
   ylabel(' Magnetic Polarity')
   title (' Magnetic Timescale for 10 Myr')
   pause
   clf
   axis;
%
%  create a timescale that is symmetric about the ridge axis
%
   [ndat,mdat]=size(t);
   mt =-fliplr(t(2:ndat)')';
   mp = fliplr(p(2:ndat)')';
   time=[mt;t];
   polarity=[mp;p];
%
% now create subset of timescale for near ridge magnetics
%
   [ndat,mdat] = size(time);
   nmid = fix(ndat/2);
%
% time is in increments of 10k, so +/- 10My will be 2048 long vector
% (FFTs are made for powers of two)
%
   ntime=2048;
   ntime2=fix(ntime/2);
   stime = time((nmid-ntime2+1):(nmid+ntime2),1);
   spolr = polarity((nmid-ntime2+1):(nmid+ntime2),1);
   plot(stime',spolr') 
   xlabel(' Million Years Before Present')
   ylabel(' Magnetic Polarity')
   axis([-10,10,-3,3])
   pause
   clf
%
%  load the Pacific Antarctic rise data.
%
   load NBP9707.xydm
   dist = NBP9707(:,3);
   magobs = NBP9707(:,4);
   plot(dist,magobs);
   xlabel('Distance (km)')
   ylabel('Magnetic Anomaly (n tesla)')
   pause
   clf
%
%  these are the model parameters that you will vary to fit the data
%  and a good starting point for the Pacific-Antarctic ridge data
%  halfrate v = 50000 m/my; dt = 20 my, L =v*dt is in meters, z = 3000 m
%  C is an arbitrary constant needed to match the amplitude
%  
   v = 47000.  ;  dt = 20;
   L = v*dt;
   z = 3000.;
   theta = -20.* pi / 180.;
   C = 3.e6;
   nx = 2048;
   dx = L/(nx-1);
   nx2 = nx/2;
%
% generate the distance vector and the wavenumbers, see Exercise 2.9
%
   x = (-nx/2:nx/2-1)*dx;
   k=((-nx/2:nx/2-1)/L)';
%
