%Model file: trapezod.m
%Comments: The model file, trapezod.m, is an example of trapezohedron
%model from Coggon (1976). This example shows how
%model parameters are to be given. The variables are Ncor D number of corners of the model; Hintn D total intensity of
%ambient magnetic induction, gamma; Hincl D inclination of
%Hintn, degrees, downward from horizontal; DeclDdeclination
%of Hintn, clockwise from north; Susc D magnetic volume susceptibility
%in units SI, a dimensionlessnumberequal to ((c)ör  -1),
%where (c)^r D magnetic permeability of the model relative to
%free space; Mstrength, Mincl, Mdecl D magnitude (gamma),
%inclination and declination (degrees), respectively, of remnant
%magnetic induction; Nf D number of faces; Fht D height of observation
%plane above origin, meters; and dens D density of
%model, g/cm3.
calgrv=1; % Change to zero if gravity field is not required
calmag=1; % Change to zero if magnetic field is not required
Ncor=26; Hintn=50000; Hincl=50; Decl=0; Susc=0.01;
Mstrength=0; Mincl=0; Mdecl=0;

%Comments: Corner is an array of x, y, z coordinates of corners
%in meters, one in each row, in a right-handed system with
%x-axis northward, y-axis eastward, and z-axis downward. Corners
%may be given in any order.
Corner = [100 0 0; 75 -75 0; 0  -100 0;  -75  -75 0;  -100 0 0;...
 -75 75 0; 0 100 0; 75 75 0; 75 0 -75; 60 -60 -60; 0 -75 -75;...
-60 -60 -60; -75 0 -75; -60 60 -60; 0 75 -75; 60 60  -60;...
0 0 -100; 75 0 75; 60  -60 60; 0  -75 75;  -60  -60 60;  -75 0 75;...
 -60 60 60; 0 75 75; 60 60 60; 0 0 100];
fht=200; dens=10;
Nf = 24;

% Add fht to depths of all corners
Corner(:,3)=Corner(:,3)+fht;

%Comments: In each row of Face, the first number is the number
%of corners forming a face; the following are row numbers of the
%Corner array with coordinates of the corners which form that
%face, seen in ccw order from outside the object. The faces may
%have any orientation and may be given in any order, but all
%faces must be included.
Face=zeros([50,9]); % Initialize a sufficiently large array
Face(1,1:5)=[4 1 2 10 9]; Face(2,1:5)=[4 2 3 11 10];
Face(3,1:5)=[4 3 4 12 11]; Face(4,1:5)=[4 4 5 13 12];
Face(5,1:5)=[4 5 6 14 13]; Face(6,1:5)=[4 6 7 15 14];
Face(7,1:5)=[4 7 8 16 15]; Face(8,1:5)=[4 8 1 9 16];
Face(9,1:5)=[4 9 10 11 17]; Face(10,1:5)=[4 11 12 13 17];
Face(11,1:5)=[4 13 14 15 17]; Face(12,1:5)=[4 15 16 9 17];
Face(13,1:5)=[4 1 18 19 2]; Face(14,1:5)=[4 2 19 20 3];
Face(15,1:5)=[4 3 20 21 4]; Face(16,1:5)=[4 4 21 22 5];
Face(17,1:5)=[4 5 22 23 6]; Face(18,1:5)=[4 6 23 24 7];
Face(19,1:5)=[4 7 24 25 8]; Face(20,1:5)=[4 8 25 18 1];
Face(21,1:5)=[4 20 19 18 26]; Face(22,1:5)=[4 22 21 20 26];
Face(23,1:5)=[4 24 23 22 26]; Face(24,1:5)=[4 18 25 24 26];

%Comments: Rectangular grid of stations for computing fields.
%Profiles are along the x-axis (north-south direction). All values
%are in meters.
s_end= -320; % Starting value of x; south end of profiles
stn_spcng = 40; % Stepsize in north direction; stn interval
n_end= 320; % Last x; maximum north coordinate
w_end= 0; % y value for westernmost profile
prof_spcng=1; % Profile spacing
e_end= 0; % y value for easternmost profile