clear
a = 6371.2;
t = 0;
r = 6371+500;
Dst = -50;
F107 = 100;

options.n_min = [ 1 14]; % n_min
options.n_max = [13 45]; % n_max
options.cord = 1;     % geocentric coordinates
options.pred = [1 1 1 1 0 0];
% complete path weree to find the file umdl.CM4
options.filename =  'q:\CM4\umdl.CM4';

% B along theta-profile
theta = [1:2:179];
phi = 180;
r = a+400; % 400 km altitude
[B_x, B_y, B_z] = cm4(t, r, theta, phi, Dst, F107, options);

subplot(3,2,1) 
source = 1; % core field
plot(theta, B_x(:,source), '-r', theta, B_y(:,source), '-g', theta, B_z(:,source), '-b')
axis([0 180 -50000 50000])
set(gca, 'Xtick', [0:30:180])
xlabel '\theta [deg]'
ylabel '[nT]'
legend('B_x', 'B_y', 'B_z')
title 'core field'

subplot(3,2,2) 
source = 2; % crustal field
plot(theta, B_x(:,source), '-r', theta, B_y(:,source), '-g', theta, B_z(:,source), '-b')
axis([0 180 -20 20])
set(gca, 'Xtick', [0:30:180])
xlabel '\theta [deg]'
ylabel '[nT]'
title 'crustal field'

subplot(3,2,3) 
source = 3; % primary magnetospheric field
plot(theta, B_x(:,source), '-r', theta, B_y(:,source), '-g', theta, B_z(:,source), '-b')
axis([0 180 -50 50])
set(gca, 'Xtick', [0:30:180])
xlabel '\theta [deg]'
ylabel '[nT]'
title 'primary magnetospheric field'

subplot(3,2,4) 
source = 4; % secondary magnetospheric field
plot(theta, B_x(:,source), '-r', theta, B_y(:,source), '-g', theta, B_z(:,source), '-b')
axis([0 180 -25 25])
set(gca, 'Xtick', [0:30:180])
xlabel '\theta [deg]'
ylabel '[nT]'
title 'secondary (induced) magnetospheric field'

subplot(3,2,5) 
source = 5; % primary ionospheric field
plot(theta, B_x(:,source), '-r', theta, B_y(:,source), '-g', theta, B_z(:,source), '-b')
axis([0 180 -50 50])
set(gca, 'Xtick', [0:30:180])
xlabel '\theta [deg]'
ylabel '[nT]'
title 'primary ionospheric field'

subplot(3,2,6) 
source = 6; % secondary ionospheric field
plot(theta, B_x(:,source), '-r', theta, B_y(:,source), '-g', theta, B_z(:,source), '-b')
axis([0 180 -25 25])
set(gca, 'Xtick', [0:30:180])
xlabel '\theta [deg]'
ylabel '[nT]'
title 'secondary (induced) ionospheric field'

set(gcf, 'PaperOrientation', 'Portrait', 'PaperType', 'A4', 'PaperUnits', ...
         'centimeters', 'PaperPosition', [0.2 0.2 20.5 29]); 
print -depsc 'CM4.eps'
