schumann_resonance.m

The Schumann resonance frequencies are a set of eigenmodes in the extremely low frequency (ELF) port

% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5 snapshot (COMSOL 3.5.0.477, $Date: 2008/08/28 16:03:15 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.5';
vrsn.ext = ' snapshot';
vrsn.major = 0;
vrsn.build = 477;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2008/08/28 16:03:15 $';
fem.version = vrsn;

% Geometry
g1=sphere3('6371000','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');
g2=sphere3('6371000+2e5','pos',{'1e5','0','0'},'axis',{'0','0','1'},'rot','0');
g3=geomcomp({g2,g1},'ns',{'g2','g1'},'sf','g2-g1','face','none','edge','all');

% Constants
fem.const = {'sig_ion','1e-8[S/m]', ...
  'R_earth','6371[km]'};

% Geometry

% Analyzed geometry
clear s
s.objs={g3};
s.name={'CO1'};
s.tags={'g3'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5, ...
                  'hmaxfac',[1,5e5,2,5e5,3,5e5,4,5e5,5,5e5,6,5e5,7,5e5,8,5e5,9,5e5,10,5e5,11,5e5,12,5e5,13,5e5,14,5e5,15,5e5,16,5e5], ...
                  'methodfac','triaf');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'ElectromagneticWaves';
appl.module = 'RF';
appl.gporder = 4;
appl.cporder = 2;
appl.sshape = 2;
appl.assignsuffix = '_rfw';
clear prop
prop.analysis='eigen';
prop.eigtype='lambda';
prop.divcond='on';
appl.prop = prop;
clear bnd
bnd.type = 'E0';
bnd.ind = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
appl.bnd = bnd;
clear equ
equ.sigma = 'sigma';
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'r','sqrt(x^2+y^2+z^2)', ...
  'sigma','sig_ion*(r-R_earth)/(r+dist1-R_earth)'};

% Descriptions
clear descr
descr.expr= {'sigma','Conductivity profile','r','Distance from Earth center'};
fem.descr = descr;

% Coupling variable elements
clear elemcpl
% Boundary distance variables
clear elem
elem.elem = 'elmapextr';
elem.g = {'1'};
clear bnd
bnd.src = {{'1'}};
bnd.ind = {{'1','2','3','4','13','14','15','16'},{'5','6','7','8','9', ...
  '10','11','12'}};
elem.geomdim = {{{},{},bnd,{}}};
elem.method = 'closest';
elem.opname = {'dist1_op'};
elem.srcmap = {'1'};
elem.dstmap = {'1'};
elem.map = {struct('type',{'unit'},'frame',{'ref'})};
elemcpl{1} = elem;
clear elem
elem.elem = 'elshapevar';
elem.g = {'1'};
clear equ
equ.var = {'dist1','sqrt((dist1_op(x)-x)^2+(dist1_op(y)-y)^2+(dist1_op(z)-z)^2)'};
equ.ind = {'1'};
elem.lagorder = '1';
elem.geomdim = {{{},{},{},equ}};
elemcpl{2} = elem;
fem.elemcpl = elemcpl;

% Descriptions
descr = fem.descr;
descr.const= {'sig_ion','Conductivity at ionosphere boundary','R_earth','Earth radius'};
fem.descr = descr;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;

% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femeig(fem, ...
               'complexfun','on', ...
               'solcomp',{'psi','tExEyEz10','tExEyEz20','tExEyEz21'}, ...
               'outcomp',{'psi','tExEyEz10','tExEyEz20','tExEyEz21'}, ...
               'blocksize','auto', ...
               'neigs',20, ...
               'linsolver','pardiso');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'slicedata',{'Wav_rfw','cont','internal','unit','J/m^3'}, ...
         'slicexspacing',5, ...
         'sliceyspacing',0, ...
         'slicezspacing',0, ...
         'slicemap','jet(1024)', ...
         'solnum',1, ...
         'title','lambda(1)=7.525985    Slice: Total energy density, time average [J/m^3]', ...
         'grid','on', ...
         'camlight','on', ...
         'campos',[-2.636014311828346E7,-3.43532073434725E7,2.4999999999999996E7], ...
         'camtarget',[0,0,0], ...
         'camup',[0,0,1], ...
         'camva',39.59775270904987);

% Plot solution
postplot(fem, ...
         'tetdata',{'normE_rfw','cont','internal','unit','V/m'}, ...
         'tetmap','jet(1024)', ...
         'tetkeep',1, ...
         'tetkeeptype','random', ...
         'solnum',1, ...
         'title','lambda(1)=7.525985    Subdomain: Electric field, norm [V/m]', ...
         'grid','on', ...
         'camlight','on', ...
         'campos',[-2.636014311828346E7,-3.43532073434725E7,2.4999999999999996E7], ...
         'camtarget',[0,0,0], ...
         'camup',[0,0,1], ...
         'camva',39.59775270904987);

% Plot solution
postplot(fem, ...
         'tetdata',{'normE_rfw','cont','internal','unit','V/m'}, ...
         'tetmap','jet(1024)', ...
         'tetkeep',1, ...
         'tetkeeptype','random', ...
         'solnum',7, ...
         'title','lambda(7)=24.435697    Subdomain: Electric field, norm [V/m]', ...
         'grid','on', ...
         'scenelight','on', ...
         'campos',[-3.63267184678768E7,-3.993975152462519E7,5.82756948608742E7], ...
         'camtarget',[0,0,0], ...
         'camup',[-0.8609995738858682,0.4567973337733771,-0.22364241463972706], ...
         'camva',11.623672480256504);