The Curlometermethod and Cluster II

Christian Maszl August 31, 2004

Contents

2

Chapter 1

Introduction

1.1 Abstract

In this project, the curlometer method was applied to FGM (Flux Gate Magnetometer) data, provided by the four CLUSTER II satellites. The primary objective was to compute the current density J⃗_{av} and ∇ · B. In order to investigate data over long periods, a \begin{matrix}c \\ ◯ \end{matrix}MatLab program was created to perform such tasks. Additionally the structure of the program was designed to allow plugins for other instruments to be implemented easily (e.g. WHISER, EDI, CIS... = multipurpose software). Finally, some basic model current sheets were tested to demonstrate the limitations of the curlometer method quite plainly.

Summary

2.1 Development of a Multipurpose Software

A multipurpose software was developed and partially tested. Unfortunately there were some restrictions, which are described in sec(4.3) on p.(28) and in sec(4.4) on p.(29). The following instruments are supported: FGM, CIS, EDI, ASPOC, EFW and WHISPER. Furthermore it is possible to obtain position- , phase- and spin axis orientation data of the spacecrafts. High resolution- and normal resolution mode for FGM magnetometer data is fully functional. The current density J⃗_{av} and ∇·B⃗ were calculated by using the Curlometer method from Sec(3.5). This can be performed for a single day or longer periods up to one year duration. The current densities J_{av} are typically in the order of 10^{−7} Am^{−2}.

2.2 Investigations on Model Current Sheets

Four main effects in case of a current tube with homogenous current density J⃗_{av} can be observed:

1) Seemingly widening of the tube diameter by increasing the satellite separation.

2) Decrease of the current density amplitude with increasing satellite separation.

3) First appearance of backward currents if s ≈ r.

4) ∇ · B⃗ not useable as an error estimate for thin current tubes.

In case of a flat current sheet, only the points 1) and 2) can be observed. Due to the special character of these sheets the next step could be, to make this investigations with more artless current sheets. It is possible that case 3) not occur in nature because of the shape of the magnetopause. Because of the large scales I estimate that the sheet is locally flat and therefore point 1) and 2) are the main sources of error.

Background

3.1 Cluster II - The Mission

Cluster^{1} is a collection of four spacecraft flying in formation around Earth. They relay the most detailed ever information about how the solar wind affects our planet in three dimensions. The solar wind (the perpetual stream of subatomic particles given out by the Sun) can damage communications satellites and power stations on Earth. The operation life-time of the Cluster mission runs from February 2001 and December 2005.

The four Cluster spacecraft have spent several years passing in and out of our planet’s magnetic field. Their mission will be to complete the most detailed investigation ever made into the ways in which the Sun and Earth interact.

3.1.1 What’s special?

The Sun emits the solar wind, which is a thin, hot, ionised gas that carries particles and magnetic fields outward from the Sun.

The Earth is shielded from the full blast by its magnetosphere, the region around our planet controlled by its magnetic field. Some solar wind descends into Earth’s upper atmosphere through the polar cusps, funnel-like openings in the magnetosphere at the poles. These energetic particles excite atoms and molecules in the upper atmosphere to create the Northern and Southern Lights (the auroras). The part of a planetary magnetosphere that is pushed in the direction of the solar wind is known as the magnetotail. Cluster will determine the physical processes involved in the interaction between the solar wind and the magnetosphere by visiting key regions like the polar cusps and the magnetotail. The four Cluster spacecraft map the plasma structures contained in these regions in three dimensions. The simultaneous four-point measurements also allow close studies of plasma quantities in both space and time. During periods of high solar activity (which cycles every 11 years), the solar wind can be particularly energetic. This can have a dramatic effect on human activities, disrupting electrical power and telecommunications or causing serious problems in the operation of satellites, especially those in geostationary orbit. Subtle changes to the weather on Earth also occur during these times. Watching the effects of this increased activity during these periods is one of the main tasks of Cluster.

1Excerpt: www.esa.int - ESA science / Cluster II - Overview

6 CHAPTER 3. BACKGROUND

Figure 3.1: The Cluster II tetrahedron Figure 3.2: Artist’s view of the solar wind and the magnetossphere

Figure 3.3: Cluster II trapping region

Understanding the interaction between the solar wind and the magnetosphere and how the plasma levels of the magnetosphere are affected is important. Cluster will help us to prepare for the effects of sudden bursts of solar energy here on Earth.

3.1.2 Spacecraft

The Cluster spacecraft resemble giant ’Lego’ sets, assembled from thousands of individual blocks. Each one is shaped like a giant disc, 1.3 metres high and 2.9 metres wide, with a cylinder in the centre. Six spherical fuel tanks are attached to the outside of this central cylinder. The fuel they carry accounts for more than half the launch weight of each spacecraft. Most of the fuel is consumed soon after launch and in complex manoeuvres to reach their operational orbits. Each spacecraft also carries eight thrusters for smaller changes of orbit. Around the central cylinder is the main equipment platform. Electrical power comes from six curved solar panels attached around the outside of the platform. Five batteries are used for power supply during the four-hour-long eclipses when the spacecraft enter Earth’s shadow. Rod-shaped booms open out once Cluster reaches orbit. There are two antennae for communications, two sensors, and four wire booms that operate when the spacecraft begins to spin. These measure changing electrical and magnetic fields around each spacecraft.

3.1.3 Journey

At each launch, two Cluster satellites were placed in an elliptical orbit whose height varied from 200 to 18 000 kilometres above Earth. The two satellites of each launch were then released, one after the other and used their own on-board propulsion systems to reach the final operational orbit (19 000 to 119 000 kilometres from the planet). The first pair of Cluster satellites lifted off on 16 July 2000, the second pair one month later. This gap allowed fewer people to be used for mission control in the European Space Operations Centre (ESOC) in Darmstadt (Germany). Once the booster reached the correct altitude, after liftoff, the Fregat payload assist module and its two Cluster spacecraft were released. The Fregat main engine fired almost immediately to achieve a circular orbit of approximately 200 kilometres high. About an hour later, the Fregat engine fired again to inject the spacecraft into an elliptical orbit. The two satellites were released, one after the other. Each Cluster spacecraft main engine performed six major manoeuvres, using the large amount of on-board fuel (about half of each satellite’s launch mass).

3.1.4 History

The Cluster mission was first proposed in November 1982. The idea was developed into a proposal to study the ’cusp’ and the ’magnetotail’ regions of the Earth’s magnetosphere with a polar orbiting mission. The Cluster idea developed into a proposal and then a mission. In 1996, Cluster was ready for launch. Cluster was expected to benefit from a ’free’ launch on the first test flight of the newly developed Ariane-5 booster. After several minor delays, Ariane-501 lifted off from Kourou, French Guiana on 4 June 1996, carrying its payload of four Cluster satellites. Unfortunately, intense aerodynamic loads resulted in its break-up and initiation of the automatic destruct system. To recover some of the unique science from the mission, ESA decided to build a fifth Cluster satellite (named ‘Phoenix’). It would be equipped with flight spares of the experiments and subsystems prepared for the Cluster mission. Phoenix was expected to be fully integrated and tested by mid-1997, opening the way for a launch later that year. However, awareness grew that the scientific objectives of the Cluster mission could not be met by a single spacecraft. There were proposals to rebuild three or four full-size Cluster spacecraft alongside Phoenix. After a preliminary study, it was decided that a Soyuz rocket could launch a pair of Cluster spacecraft. However, the very eccentric orbit required a new upper stage. Two flights were successfully done at the beginning of 2000 and about 6 months later, Cluster was launched by a Soyuz-Fregat launcher from Baikonur Cosmodrome, Kazakhstan.

3.1.5 Partnerships

Prime contractor for the original (lost) Cluster and replacement Cluster satellites was Dornier Satellitensysteme GmbH (now Astrium), Friedrichshafen, Germany, the leader of an industrial consortium involving 35 major contractors from all of the ESA member countries and the United States. Each spacecraft carries an identical set of 11 instruments to investigate charged particles, electrical, and magnetic fields. These were built by European and American instrument teams led by Principal Investigators. The Cluster scientific community includes the ESA Project Scientist, 11 Principal Investigators, and more than 250 Co-Investigators from ESA Member States, the United States, Canada, China, the Czech Republic, Hungary, India, Israel, Japan, and Russia.

3.2 Cluster II - Operation & Instruments

The four Cluster^{2} spacecraft are placed in nearly identical, highly eccentric polar orbits. In the two-year working period the satellites are able to make a detailed exploration of all interesting areas of the magnetosphere.

The orbit for each spacecraft is selected so that each is located at a vertex of a pre-determined tetrahedron when crossing the regions of interest within the magnetosphere. The separation distances between the spacecraft will be adjusted during the mission depending on what will be studied and will vary from a few hundred kilometres to a few Earth radii.

In order to perform measurements the spacecraft are gyrating around their z-axis. In normal operation mode the frequency is 0.25Hz in high resolution ”burst” mode the gyration frequency is 22Hz.

3.2.1 FGM - Fluxgate Magnetometer

Consists of two magnetometers^{3} on a five-metre long boom from the spacecraft to avoid interference. It measures the magnetic field along the orbit. The fluxgate magnetometers are similar to many previous instruments flown in Earthorbit and on other planetary and interplanetary missions. In order to minimise the magnetic background of the spacecraft, one of the magnetometer sensors (the outboard, or OB sensor) is located at the end of one of the two 5.2m radial booms of the spacecraft, the other (the inboard, or IB sensor) at 1.5m inboard from the end of the boom. The boom mounted sensor units are essentially passive devices, consisting of coils of wire wrapped around toroidal cores of magnetic material. The magnetometers have eight possible operating ranges; of these, five are used on the Cluster magnetometers. These ranges were selected to provide good resolution in the solar wind (with expected field magnitudes between 3 and 30nT), and up to the highest field values expected in the magnetosphere along the Cluster orbit (up to about 1000nT). The highest range (≈ 65500nT) is used only to facilitate background testing.

3.2.2 EDI - Electron Drift Instrument

The Electron Drift Instrument (EDI) measures the drift of a weak beam of test electrons, from which the ambient electric fields and magnetic field gradients can be determined. As a by-product, the magnetic field strength is also measured.

3.2.3 ASPOC - Active Spacecraft Potential Control Experiment

Neutralises the positive charge that is built on the spacecraft which would seriously interefere with the measurements. It emits indium ions to space through a small needle.

2Excerpt: spaceweb.oulu.fi/projects/cluster/ - Cluster II project in Oulu

3Excerpt: www.sp.ph.ic.ac.uk/Cluster/ - Imperal College Cluster Group

10 CHAPTER 3. BACKGROUND Figure 3.4: Anatomy of the Cluster Satellites

3.2.4 STAFF - Spatio-Temporal Analysis of Field Fluctuation experiment

Consists of a magnetometer on the end of a five-meter long boom that analyzes waves (rapid variations in the magnetic field).

3.2.5 EFW - Electric Field and Wave Experiment

Measures electric field to study plasma convection and waves with sensors on four 50 m long wire booms.

3.2.6 DWP - Digital Wave Processing experiment

Controls the wave experiments and performs calculations for them. Helps the WECinstruments to make effective use of spacecraft power resources and telemetry information.

3.2.7 WHISPER - Waves of High frequency and Sounder for Probing of Electron density by Relaxation experiment

Measures the density of charged particles that fill near-Earth space with resonance sounding (radar) technique.

3.2.8 WBD - Wide Band Data instrument

Makes high-resolution measurements of both electric and magnetic fields in selected frequency bands from radio sounds from particles around Earth’s magnetic poles.

3.2.9 PEACE - Plasma Electron And Current Experiment

Looks at all electrons in space which have low to medium energies, counts them and measures their direction and speed. One part of the instrument measures faster electrons, another part slower ones.

3.2.10 CIS - Cluster Ion Spectrometry experiment

Analyses the composition, mass and distribution functions of ions in the nearby space plasma and in the solar wind during each four-second spin of the spacecraft.

3.2.11 RAPID - Research with Adaptive Particle Imaging Detectors

Uses two different and independent systems (Imaging Electron Spectrometer and Imaging Mass Spectrometer) to record the highest energy electrons and ions.

3.2.12 WEC - Wave Experiment Consortium

DWP, EFW, STAFF, WBD, and WHISPER, which all measure electric and magnetic fields and waves in some form. They have been grouped together to form WEC in order to maximise resources

Figure 3.5: Topography of the solar-terrestrial enviroment

3.3 The Physical Background

3.3.1 Solar Wind

The sun^{4} emits a highly conducting plasma at supersonic speeds of about 500km/s into the interplanetary space. The solar wind consits mainly of electrons and protons, with an admixture of 5% Helium. Because of the high conductivity, the solar magnetic field is frozen in the plasma and drawn outward by the expanding solar wind. Typical values of solar wind parameters: n_{e} ≈ 5cm^{−3}, T_{e} ≈ 10^{5}K (1eV = 11600K) and B of the interplanetary magnetic field is of the order of 5nT.

When the solar wind hits on the Earth’s dipolar magnetic field, it cannot simply pentrate it bur rather is slowed down and, to a large extent, deflected around it. Since the solar wind hits the obstacle with supersonic speed, a bow shock wave is generated, where the plasma is slowed down and a substanial fraction of the particels’ kinetic energy is converted into thermal energy. The region of thermalized subsonic plasma behind the bow shock is called magnetosheath. Its plasma is denser and hotter than the solar wind plasma and the magnetic field strength has higher values in this region.

3.3.2 Magnetossphere

The interplanetary magnetic field lines cannot pentrate the terrestial field lines and the solar wind particles cannot leave the interplanetary field lines due to the frozen-in characteristic of a highly conducting plasma.

4Excerpt: W. Baumjohann / R. A. Treumann - Basic Space Plasma Physics

Figure 3.6: Plasma structure of Earth’s magnetosphere

The boundary separating the two different regions is called magnetopause and the cavity generated by the terrestrial field has been named magnetosphere. The kinetic pressure of the solar wind plasma distorts the outer part of the terrestrial dipolar field. At the frontside the field becomes compressed, while the nightside magnetic field is stretched out into a long magnetotail which reaches far beyond the lunar orbit.

3.3.3 Radiation belt

The radiation belt lies on dipolar field lines between about 2 and 6R_{E} (1R_{E} ≈ 6371km). It consists of energetic electrons and ions which move along the field lines and oscillate back and forth between the two hemispheres. This is due to the fact that the dipolar magnetic field acts like a magnetic mirror. A dipole magnetic field has a field strength minimum at the equator and converging field lines at the north and south pole. Since the magnetic moment is invariant and the total energy is a constant of motion, only the pitch angle can change when the magnetic field increases or decreases along the guiding center trajectory.

μ = \frac{mv^{2}sin^{2}α}{2B}  v_{⊥} = vsinα  (3.1)

Eq(3.1) shows that the pitch angles of a particle at different locations are directly related to the magnetic field strengths at those locations according to:

\frac{sin^{2}α_{2}}{sin^{2}α_{1}} = \frac{B_{2}}{B_{1}}  (3.2)

If the pitch angle reaches α = 90, the particle is reflected from this mirror point. The bounce period τ_{b} is the time it takes a particle to move from the equatorial plane to one mirror point, then to the other and back to the equatorial plane.

τ_{b} = 4 ∫_{0}^{λ_{m}} \frac{ds}{v_{∥}}  (3.3)

The bounce period of 1keV electrons is in the order of a few seconds while it takes a few minutes for the heavier protons to complete a full bounce cycle.

Typical parameters: n_{e} ≈ 1cm^{−3}, T_{e} ≈ 5 · 10^{7}K and B ≈ 100 − 1000nT

3.3.4 Plasma sheet

Most of the magnetotail plasma is concentrated around the tail midplane in an about 10R_{E} thick plasma sheet. Near the Earth, it reaches down to the high-latitude auroral ionosphere along the field lines. The plasma sheet is also the source of the ring current.

Typical parameters: n_{e} ≈ 0.5cm^{−3}, T_{e} ≈ 5 · 10^{6}K and B ≈ 10nT

3.3.5 Magnetotail lobe

The magnetotail lobe contains a highly rarefied plasma.

Typical parameters: n_{e} ≈ 10^{−2}cm^{−3}, T_{e} ≈ 5 · 10^{5}K and B ≈ 30nT

3.3.6 Ionosphere

The solar ultraviolet light impinging on the Earth’s atmosphere ionizes a fraction of the neutral atmosphere. At altitudes above 80km, collisions are too infrequent to result in rapid recombination and so a permanent ionized population called the ionosphere is formed.

Typical parameters: n_{e} ≈ 10^{5}cm^{−3}, T_{e} ≈ 5 · 10^{3}K and B ≈ 10^{4}nT

The ionosphere extends to rather high altitudes and, at low- and mid-latitudes, gradually merges into the plasmasphere. The plasmasphere is a torusshaped volume inside the radiation belt. It contains a cool but dense plasma of ionospheric origin, which co-rotates with the Earth.

Typical parameters: n_{e} ≈ 5 · 10^{2}cm^{−3}, T_{e} ≈ 5 · 10^{3}K

The plasmasphere extends outwards to about 4R_{E} where the density drops down sharply to about 1cm^{−3}. This boundary is known as the plasmapause. At high latitudes plasma sheet electrons can precipitate along magnetic field lines down to ionospheric altitudes where they collide with and ionize neutral atmospheric particles. As a by-product, photons emitted by this process create the polar light known as the aurora. These auroras are typically observed inside the auroral oval which contains the footprints of those field lines which thread the plasma sheet. Inside of the auroral oval lies the polar cap which is threaded by field lines connecting to the tail lobe.

Figure 3.7: Synopsis of magnetospheric currents

3.4 Magnetospheric currents

The plasmas discussed in the previous section are usually not stationary but move under the influence of external forces. Sometimes ions and electrons move together, like in the solar wind. But there are other occasions and in other plasma regions where ions and electrons move in different directions creating an electric current. Such currents are very important for the dynamics of the Earth’s plasma environment. They transport charge, mass, momentum and energy. Moreover, these currents create magnetic fields which may severely alter or distort any pre-existing fields.

In fact, the distortion of the terrestial dipole field into the typical shape of the magnetosphere is accompanied by electrical currents. The compression of the terrestrial magnetic field on the dayside is associated with current flow across the magnetopause surface, the magnetopause current. The taillike field of the nightside magnetosphere is accompanied by the tail current which flows on the tail surface and the neutral sheet current in the central plasma sheet. Both of these currents are connected and form a θ-like current system if viewed from along the Earth-Sun line.

Another large-scale current system, which influences the configuration of the inner magnetosphere is the ring current. The ring current flows around the Earth in a westward direction at radial distances in the order of several Earth radii and is carried by the radiation belt particles. In addition to their bounce motion these particles drift (due to the curvature of the magnetic field) slowly around the Earth. Since the protons drift westward while the electrons move in the eastward direction,

Figure 3.8: Specular reflection off the magnetopause

this constitutes a net charge transport.

A number of current systems exist in the conducting layers of the Earth’s ionosphere at altitudes of 100-150km. Most notable are the auroral electrojets located inside the auroral oval, the Sq currents in the dayside mid-latitude ionosphere, and the equatorial electrojet near the magnetic equator.

In addition to these perpendicular currents, currents also flow along magnetic field lines. The field-aligned currents connect the magnetospheric current systems in the magnetosphere to those flowing in the polar ionossphere. The field-aligned currents are mainly carried by electrons and are essential for the exchange of energy and momentum between these regions.

3.4.1 Magnetopause current

The following chapter deals with the magnetopause current. The fully ionized and magnetized solar wind plasma cannot mix with the terrestrial flux tubes. Instead, it will deviate from its original direction and will, by its dynamical pressure, compress the terrestrial field and confine it into a small region of space, producing the magnetosphere. This layer is a discontinuity which must be different from bow shock because the plasma flow behind the bow shock is subsonic. In fact, to first order, the magnetopause can be regarded as a tangential discontinuity. Separating the solar wind from the magnetospheric magnetic field and being a surface across which the magnetic field strength jumps from low interplanetary value to high magnetospheric field strength, the magnetopause represents a surface current

Figure 3.9: 3D geometry of magnetopause currents

layer. The origin of this current is shown schematically in fig(3.8). Specularly reflected ions and electrons hitting the magnetospheric field inside the magnetopause boundary will perform half a gyro-orbit inside the magnetic field before escaping with reversed normal velocity from the magnetopause back into the magnetosheath. The thickness of the solar wind-magnetosphere transition layer under such idealized conditions becomes of the order of the ion gyroradius.

r_{gi} = \frac{v_{sw⊥}}{ω_{gi}}  (3.4)

Where v_{sw⊥} is the normal velocity component of the solar wind and ω_{gi} represents the ion cyclotron frequency given by:

ω_{gi} = \frac{q_{i^{+}}B}{m_{i^{+}}}  (3.5)

For typical values of B = 10nT and v_{sw⊥} = 300m/s one gets a value of r_{gi} = 1000km for the thickness. Electrons also perform half gyro-orbits, but with much smaller gyroradii. The sense of gyration inside the boundary is opposite due to the sign of the charge for both kinds of particles leading to the generation of a narrow surface current layer. This current provides the additional magnetic field which compresses the magnetospheric field into the magnetosphere and at the same time annihilates its external part. It is a diamagnetic current caused by the perpendicular density gradient at the magnetopause. The current density inside the magnetopause can be estimated to approximately 10^{−6} Am^{−2}. The total current flowing inside the magnetopause is of the order of 10^{7}A. In the equatorial plane the magnetopause current flows from dawn to dusk, as shown schematically in fig(3.9). It closes on the tail magnetopause, where it splits into northern and southern parts flowing across the lobe magnetopause from dusk to dawn. The tail magnetopause current is additionaly fed by the cross-tail neutral sheet current which flows from dawn to dusk.

Figure 3.10: Geometrical coherences of Cluster II

3.5 The Curlomether Technique

3.5.1 Derivation of J⃗_{av} and ∇ · B⃗

Assuming that the current density J⃗_{av} is constant over the tetrahedal volume and the magnetic field B⃗ changes only linearly, one can use Ampere’s law to compute J⃗_{av}. The difference from ∇ · B⃗ = 0 gives a quality estimate of the derived quantities. Geometrical coherences can be observed in fig(3.10).

μ_{0}J⃗_{av} = ∇ × B⃗  (3.6)

Applying the theorem of Stoke’s and assuming that J⃗_{av} is constant over the tetrahedral volume one gets

μ_{0}J⃗_{av} ∫_{A} d⃗a = ∫_{L} d⃗lB⃗  (3.7)

Assuming that B⃗ changes only linearly and implying that the satellite separation is small compared with the magnitude of the variation of B⃗ one can write (using SC1 as reference)

B⃗_{12} = \frac{B⃗_{1} + B⃗_{2}}{2}  r⃗_{12} = r⃗_{1} − r⃗_{2}  (3.8)
B⃗_{23} = \frac{B⃗_{2} + B⃗_{3}}{2}  r⃗_{23} = r⃗_{2} − r⃗_{3}  (3.9)
B⃗_{13} = \frac{B⃗_{1} + B⃗_{3}}{2}  r⃗_{13} = r⃗_{1} − r⃗_{3}  (3.10)

Rewritting

∫_{A} d⃗a  as  \frac{1}{2}(r⃗_{13} × r⃗_{12})

and projecting the magnetic field B⃗ on the edges of the tetrahedron one gets

\frac{1}{2}μ_{0}J⃗_{av} · (r⃗_{13} × r⃗_{12}) = B⃗_{12} · r⃗_{12} − B⃗_{23} · r⃗_{23} − B⃗_{13} · r⃗_{13}.  (3.11)

After some mathematical transformations using eq(3.8) - eq(3.10) and defining ΔB⃗_{13} = B⃗_{1} − B⃗_{3}, ΔB⃗_{12} = B⃗_{1} − B⃗_{2} one finally achieves

μ_{0}J⃗_{av} · (r⃗_{13} × r⃗_{12}) = ΔB⃗_{13} · r⃗_{12} − ΔB⃗_{12} · r⃗_{13}.  (3.12)

In order to obtain a quality estimate of the computed current density J⃗_{av}, one investigates the aberration of ∇ · B⃗ from 0. J⃗_{av} and ∇ · B⃗ are subject to errors of which there are basically three types. The first concerns to the measurement uncertainties in B⃗ and the spatial configuration of the four spacecrafts.The second relates to the linear interpolation which is made between the various measurement points. The third relates to the simultaneity of the measurements.

Applying the theorem of Gauss on ∇ · B⃗ one gets

∫_{V} dV ∇ · B⃗ = ∫_{A} d⃗a · B⃗.  (3.13)

Treating B⃗ as constant and taking the average of ∇ · B⃗ yields

< ∇ · B⃗ >_{av} V = | ∑_{ν=1}^{4} A⃗_{ν} · B⃗_{mν} | .  (3.14)

ν represent the four areas of the tetrahedron. For example A⃗_{1} one gets

A⃗_{1} = \frac{1}{2}(r⃗_{13} × r⃗_{12}).  (3.15)

For the tetrahedral volume one can write (using SC1 again as reference)

V =| \frac{1}{3!} r⃗_{14} · (r⃗_{13} × r⃗_{12}) | .  (3.16)

At this state we use the linear approximation for B⃗. For the triangel composed of spacecrafts 1, 2 and 3 one can, for the interception point of the three heights, write

B⃗_{m1} = \frac{B⃗_{1} + B⃗_{2} + B⃗_{3}}{3} .  (3.17)

Analogous to this, one can obtain the fields for the remaining three areas. Inserting eq(3.15), eq(3.16) and eq(3.17) into eq(3.14) one obtains the result (α, β, γ = 2, 3 and

< ∇ · B⃗ >_{av} r⃗_{14} · (r⃗_{13} × r⃗_{12}) = | ∑_{cyclic} ΔB⃗_{1α} · (r⃗_{1β} × r⃗_{1γ}) |.  (3.18)

To get an error estimate for the computed current density one takes eq(3.18) divided by μ_{0} and divides this by the norm of eq(3.12). One performs the first division to achieve the same dimension. Finally one can write for the error estimate

{P  = \atop μ_{0} \frac{< ∇ · B⃗ >_{av}}{| J⃗_{av} |}.}  (3.19)

3.5.2 Derivation of J_{av⊥} and J_{av∥}

In order to investigate the currents created by ions and electrons we split J⃗_{av} into perpendicular and parallel parts with respect to the magnetic field strength B⃗. The field aligned currents J_{av∥} are mainly carried by electrons and are essential for the exchange of energy and momentum between the magnetospheric current systems in the magnetosphere and in the polar ionosphere. The perpendicular currents J_{av∥} are carried by the ions. First we calculate the angle θ between J⃗_{av} and B⃗,

cos θ = \frac{B⃗ · J⃗_{av}}{| B⃗ || J⃗_{av} |}.  (3.20)

Splitting J⃗_{av} into a perpendicular and parallel part one derives

J_{av∥} =| J⃗_{av} | cos θ  J_{av⊥} =| J⃗_{av} | sin θ.  (3.21)

Using the fact that

sin θ = \sqrt{1 − cos^{2} θ}  (3.22)

and taking into account that \begin{matrix}c \\ ◯ \end{matrix}MatLab is more rapid with matrices, one substitutes cos θ and finally achieves for J_{av∥}

J_{av∥} = \frac{B⃗ · J⃗_{av}}{| B⃗ |}.  (3.23)

Similarly one derives for the perpendicular part J_{av⊥}

J_{av⊥} =| J⃗_{av} | \sqrt{1 − (\frac{B⃗ · J⃗_{av}}{| B⃗ || J⃗_{av} |})^{2}}.  (3.24)

The Software

4.1 Preface

Originally this program was planed to apply the curlometer method to FGM magnetometer data for long periods. The first version was capable of searching for ascii data files and performing the calculation. In order to save processing time and improve the flexibility of the program, the more flexible subroutine getData of Yuri Kothyaintsev (IRFU) was implemented as raw data source. A disadvantage incurred as a result of this step disallows the program to be run locally on a PC. As a consequence, it is necessary to establish a SSH connection to SANNA and run the program there. Due to the high flexibility of the routine getData and the ability to obtain data for many measurement instruments (EFW, EDI, WHISPER ...) Stefan Buchert (IRFU) suggested that the program be extended. This would facilitrate other members of the institute which could write plugins for their special tasks. As a result the program was altered so that it now searches, obtains and sorts data for later disposal. Furthermore it provides a list with errors, gaps, ... that might occur during the operation.

Because of the programs ”history” it is possible to find structures in the code that are not necessary anymore. e.g. The first approach to fetch data from binary file was to use several c-programs of Stefan Buchert (IRFU) and the use of pipes to get the data. Therefore the subroutien filelist search for binary raw data files in automatic mode. After the switch to getData this isn’t necessary anymore because this routine need only the date as input parameter. GetData fetches the data from ISDAT database CDF files.

It is possible to find more of such parts. These parts should be removed in a maybe second revision. In the following section you will find information about all used subroutines, limitations, bugs and not finished tasks. There are also a lot of coments in the source code to ease debugging and to reduce the period of vocational adjustment.

4.2 Description of Program Parts

4.2.1 The Main Program - mainprogram.m

This section provides a rough overview of the main program. For details relating to each subroutine see the following chapters.

Figure 4.1: Artist’s view of Cluster II and Earth

syscheck.m

The entire program is controlled by a small configuration file configfile.m. Initially the subroutine syscheck.m attempts to check the entries of configfile.m, next it performs several write and read tests and finally sets variables for the different operation modes.

filelist.m

The next step is the routine filelist.m. This searches for datafiles in automatic mode. In manual mode this routine is skipped. Due to the fact that getData doesn’t require a filelist anymore this subroutine can be replaced with a program that generates a date list for a whole year. If time is short or there is a lack of ”motivation” to program it, it is possible to set the variable mode in syscheck.m to fn. This will generate a list of FGM raw data files. getData will take the date out of the filenames and so will work properly.

userinput.m

In automatic mode userinput.m asks for the desired operation mode. The first possibility is choice [A]. Following selection, hte program performs calculation for an entire year of data. This can require approximately 4h for normal mode-4s data or a few days for high resolution-45ms data depending on the processor speed and the available RAM (The datafiles are quite large and consume much memory. Therefore, it would be recommended to use a system with at least 512MB to avoid ”swapping”.)

The choice item [F] gives a filelist of found data and [M] allows to define a manual range. In manual mode the desired date can be selected.

yuri.m

This subroutine provides the raw data and information concerning the ”fetching” process. This subroutine is the interface to getData.m.

sorting.m

The four satellites do not function precisely in synchrony. Therefore sorting.m synchronises the four data files together depending on the time resolution of the input files.

report.m

report.m stores a report file to the HD in case yuri.m relate some errors.

jav fgm.m

jav fgm.m applies the curlometer method to the data and saves the raw data and the computed data to the HD.

4.2.2 The Control Center - configfile.m

is the main control unit of the program. In order to reduce the number of dialogues during execution of the program all necessary parameters become stored in this file.

datapath

This variable stores the path to the input data files.

Default on SANNA: ’/data/cluster/DDS’

isdat db

ISDAT database string.

Default on SANNA: ’disco:10’

year

defines the year to analyse for the automatic mode. Manual mode overrides this with a user entry.

op mode

selects the instrument. At the moment only ’nr’ and ’hr’ modes (FGM - 4s and 45ms data) are tested and working. extension

is a relict from a early development stage. Move to syscheck.m to obviate user errors if you like.

Default: ’*’

datatype

defines to kind of the stored datafiles. It’s possible to choose between binary ’-mat’ files and with any editor readable ’-ascii’ files. The recomendation is clear on binary files because loading and saving is about five times faster than in ascii format.

Default: ’-mat’

logfilepath

is the directory were the error.log and myfilelist.log created by logfiles.m in filelist.m are going to be stored. Make sure that you have write permission for this directory!

outputspath

is the directory were the computed data and the raw data is going to be stored. Make sure that you have write permission for this directory!

tempdirpath

is the directory for temporary data created by getData.m. Make sure that you have write permission for this directory!

y0

is the magnetic permeability constant μ_{0} in SI units.

Default: μ_{0} = 1.256637061 · 10^{−6} \frac{V s}{Am}

format long

is the way how \begin{matrix}c \\ ◯ \end{matrix}MatLab displays data on the screen.

advanced

The program is intended to be used for whole days in manual mode or longer periods in automatic mode. Altough it is possible to analyse smaller periods. This handles the variable advanced. The first three cells of the array represent the starttime of the computation in hours, minutes and seconds. The fourth cell is the duration. Make sure that the period isn’t longer than the remaining day.

Default: [00 00 00 86400]

compmode

swiches between manual mode (’m’) and automatic mode (’a’). nolog

set to ’n’ to surpress the generation of logfiles.

Default: ’y’

plot mode normal - Graphics output

is one of the control variables for the graphical output. Note that plotdiagram.m work only for the modes ’nr’ and ’hr’! Set this variable to ’y’ to display Jav_{x}, Jav_{y}, Jav_{z} and \frac{∇·B⃗}{|∇×B⃗|}. Choose ’n’ to skip this plot.

Default: ’y’

plot mode parnor - Graphics output

set to ’y’ produce the plot with Jav_{⊥}, Jav_{∥} and \frac{∇·B⃗}{|∇×B⃗|}.

Default: ’y’

legend on - Graphics output

toggles the legend in the graph B_{x}, B_{y}, B_{z} and | B⃗ | on or off.

Default: ’y’

4.2.3 The Security Unit - syscheck.m

This program loads configfile.m and check if the values in the configfile are valid and perform several write tests . It also converts the mode variable of the configfile to the proper input parameters for the subroutine getData.m.

At the beginning of the program one find several error messages that will be displayed in an error case. This part is followed by

Conversion to satellite instrument

This section contains the variables which are of no interest for the user. Before inserting a new plugin (e.g for EDI) make sure that the variables are checked and fit to your request.

Raw data filenames are usually in the following shape : [date(yymmdd)][change number][mode].[satellite] - e.g 0203040000fn.01. This is a magnetometer raw datafile for 4^{th} of March 2002 for satellite nr.1 with no changes. The variable mode is used in filelist.m to know for which data to look (e.g FGM - mode=’fn’).

instrument is the parameter for getData.m which selects the operation mode. You find a list of all possible parameters at the end of mainprogram.m or yuri.m. There you also find the style of the ”fetched” data.

option is also a parameter for getData.m and suppres the storage of the fetched data to the HD. The default option is ’nosave’ because the data will be stored afterwards anyway. (The only execption is the ’ve’-mode. Check syscheck.m for the details.

timeacc is the accuracy parameter for sorting.m. This parameter should be half the time resolution of the data. Note that a to low accuracy will lead to wrong results in computation and a to high accuracy will eliminate a lot of data.

timeres, datasrc, dataname and varname are used to display messages on the screen.

Checks and Tests

The rest of the subroutine is dedicated to checks, tests and also makes sure that all path details end with a /.

4.2.4 Search for Data - filelist.m

As explained above this subroutine isn’t absolute necessary any more. But it can be usefull when you are going to substitute getData.m.

In automatic mode filelist.m searches for datafiles defined in the directory datapath with the adequate ”string” defined in the mode variable. The data becomes checked and sorted. Later the programm checks if there are four files (satellite 1 to 4) for each date. All valid input files become stored in myfilelist all others in errorlog. The first part of errorlog contains files with nonmatching mode criteria. The second contains information of uncompleted quad sets. filelist.m also stores errorlog and myfilelist to the harddisk (datatodisk) in the directory defined in logfilepath (error.log and myfilelist.log).

4.2.5 Only one Question - userinput.m

In manual mode the user have to input the desired date in automatic mode the user can choose between [A] for a whole year, [M] for a favoured range or [F] to display a list of all found data files.

4.2.6 The Heart - yuri.m

This is the interface to the program getData.m. Unfortunately the subroutine get-Data.m isn’t working absolutely reliable. Check Sec(4.3) and Sec(4.4) for known bugs and recomended workarounds.

Referring to the conversation with Yuri Kothyainsev (IRFU) I recommend to substitute getData.m by a new program!

getData.m is supposed to used only for minute intervals and not hours. The main problem is that getData.m creates ascii files and pipes the data through several programs. This consume a lot of memory and time, especially when one does this for a whole day.

Prepare the input parameters

First the program take the date out of the filenames, start-times and duration out of the variable advanced. Due to the fact that getData.m can’t handle long time intervals, the duration becomes split into four pieces and four new starttimes are going to be calculated.

Next the operation mode is going to be determined. getData.m returns fetched data in cell arrays. This data is different in the number of double arrays. Therefore it’s necessary to split yuri.m in three parts. Part one for ”one” array data, part two for ”two” arrays data and part three for ”seven” arrays data, all stored in one cell.

Check the end of yuri.m or mainprogram.m for a list of all supported data files. One also finds there additional information about the cells, variable names, shape and dimension of the arrays, ...

Cells with one enclosed matrix

This is the only mode which is debugged and should work properly. Anyway there is a problem in automatic mode. If one performs calculations for longer periods then thirty days, getData.m reports ’No Data from Server’ when fetching position data. See Sec(??) for a possible workaround. First getData fetch B⃗ field data and r⃗ positon data for all satellites. It is possible that there are gaps for several hours in the data. In this case B⃗ or r⃗ are set to zeros. Success or failure of this process is reported in data B and data r. If there is no data for one satellite over a whole day, one find this information in b report and c report. a report, d report and e report contain information if there is B⃗ data and no r⃗ ...

Interpolate the position vector - av interp.m

The time resolution of the position vector r⃗ is worse than the vector B⃗. Therefore one has to interpolate r⃗ to the timeline of B⃗. This task can be done with av interp.m. The default method in yuri.m is spline.

Make clean and convert time - fromepoch.m

getData.m stores data during the process in a temporary directory defined by tempdirpath. These files are going to be cleared after usage. getData.m provides the time in ISDAT epoch. With the function fromepoch.m this format is converted into hours of a day. getData.m has also the habit to mix data sometimes if fetching data for long periods. To be sure that the timeline is in the proper order, the data becomes sorted.

Cells with two enclosed matrices

The mode of operation is similar to that with only one enlosed matrix. Check Sec(4.3.3) on p.(28) for the differences.

Cells with seven enclosed matrices

The mode of operation is similar to that with only one enlosed matrix. Check Sec(4.3.4) on p.(29) for the differences and unfinished features.

4.2.7 Synchronise the Satellites - sorting.m

This script sort the data of the four satellites according to the time. timeacc set the maximum time difference between the measurements of one event. The procedure is as following. Because of different start times of each satellite, sorting.m looks for the latest start entry. Afterwards it takes this value and compare it with the other three data files. In case a match inside the ±timeacc borders is found, the entries become stored in a list.

4.2.8 The Sneak - report.m

uses the report variables from yuri.m to generate error messages and to store them into a file. The filename is in the following shape [yymmss]0000[mode].report (e.g 0203040000fn.report). This file is going to be loaded in case of execution of plotdiagram.m.

I recommend to display this message file also for other possible plugins!

It might happen that getData.m don’t provide position data r⃗. But curlB.m will perform calculations and produce output anyhow if there is r⃗ data or not.

4.2.9 The Working Horse - jav fgm.m

contains the program curlB.m and datatodisk.m. curlB.m apply the curlometer method to the FGM data. datatodisk stores this data and the raw data to the HD. The filenames are in the following style. For the raw data [yymmdd]0000[mode].D[1-4].gse.hr.igm[asc or mat], for the current density [yymmdd]0000[mode].jav.gse.hr.igm[asc or mat] and for ∇ · B⃗ [yymmdd]0000[mode].divB.gse.hr.igm[asc or mat]. The data becomes stored in the directory defined by outputspath in configfile.m.

4.2.10 The Graphical Output - plotdiagram.m

plots diagrams on the screen. Note that this program only work for the FGM modes nr and hr. Again this program starts with syscheck.m. Therefore make sure that the input parameters in configfile.m are adjusted to your needs! plotdiagram.m is designed to plot three types of diagrams. The three diagram types first B_{x}, B_{y}, B_{z} and | B⃗ |, second Jav_{x}, Jav_{y}, Jav_{z} and \frac{∇·B⃗}{|∇×B⃗|} and third Jav_{⊥}, Jav_{∥} and

\frac{∇·B⃗}{|∇×B⃗|}. Later can be selected by changing plot mode parnor and plot mode normal in configfile.m

4.3 Unfinished Features

4.3.1 No time line for op mode=sa

Take the time lines from the position vectors and use them as first columns (yuri.m).

4.3.2 Automatic accuracy for sorting.m

The time accuracy manually defined in syscheck.m in the variable timeacc should be found by the program itself. This can be done after getData.m when the timesteps can be observed. Note that the time accuracy should be half the time resolution.

4.3.3 Test yuri.m for cells with two enclosed matrizes.

This mode isn’t fully tested yet. In order to store two arrays in one, both have to have the same length. It is possible that one of these two is zero. Therefore it is necessary to fill the empty array with the timeline of the full array and zeros for the data. This mode is tested for the EFW instrument (op mode=’ef’ in configfile.m) in manual mode for the day 030702 (yymmdd) and working. Note that we1p34 for SC1 and we4p34 for SC4 is 0.

4.3.4 Finish yuri.m for cells with seven enclosed matrizes.

This part isn’t finished yet. The only thing that is missing is to find the empty variables and to fill them with a timeline and zeros. I recommend to insert this function at line 473

4.3.5 Generate a date list for yuri.m

Due to the fact that myfilelist.m is not necessary anymore this time consuming routine can be skipped. All other routines only take the date out of the filenames. Therefore a data list should be created which can fullfill this task in a faster way. Store the date list in myfilelist.

4.3.6 Easier control of userinput.m

If one likes to enter a manual range one has to input numbers connected to the date and therefore always to check the filelist. A better way is to allow the user to define the range by declaring start- and end-times. In case myfilelist will be subsituted by a date list, the option [F] filelist can be removed.

4.3.7 Multi coordinate system support for datatodisk.m

The default filenames for datatodisk contain gse for the coordinate system and hr for the time base. The program getData.m is able to fetch data for different coordinate systems. Therefore I recommend to add information about the time unit and coordinate system in syscheck.m to have an easy way to adjust the filenames.

Note that you also have to change plotdiagram.m in this case!

4.3.8 Automatic eps-file generator for plotdiagram.m

It might be convenient, if you observe a lot of data, to create eps-files of the graphs on the fly for later use. This should avoid long loading and waiting times if you watch the data more often. This feature could also be controlled with configfile.m.

4.3.9 Replace getData.m

Due to the fact the getData.m isn’t reliable and not intended for use in the required way, I recommend to write your own data fetching routine. For FGM data you can use the C programms of Stephan Buchert (IRFU) for example. But there are also signs that the errors and problems occur in the ISDAT database, Sec(??) and Sec(??). If you have the skills you can try and find it there.

4.4 Known bugs

In this section you find a list of known bugs. Some of them are corrected. Anyhow they are listed to prevent you from problems in case you change or add program parts.