CENTRE DE RECHERCHES EN PHYSIQUE DES PLASMAS (CRPP) Association EURATOM - Confédération suisse

TCV Diagnostics

Documentation

Fast-Algorithm Bolometric Computer Aided Tomography (FABCAT)

Christian Schlatter^{1}, Jan Mlynář^{2}

CRPP EPFL PPB CH-1015 Lausanne

February, 2002

Made with L^{A}TEX

1TP IV in physics, christian.schlatter@epfl.ch 2jan.mlynar@epfl.ch

Contents

1 Introduction 1

2 Bolometry on Tokamak à configuration variable (TCV) 2

3 Organization of the code 3 3.1 Main program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.2 Supporting routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3.2.1 fabsetchord.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 3.2.2 fabsetetendue.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 3.2.3 fabsett.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 3.2.4 Bolometry movie player . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 3.2.5 Other tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

4 Operating instructions 6

4.1 How to execute

’fabcat.m’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 4.2 How to configure

’fabcat.m’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 4.3 How to read the

MDS tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

5 Detailed description of the method 10 5.1 fabcat.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 5.2 fabsig.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

5.2.1 The brightness f_{ℓ} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.2.2 Reading the MDS tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.2.3 LIUQE times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.2.4 NORMAL times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.2.5 Automatic channel

deleting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 5.2.6 Signal processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

5.3 catrec.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

5.3.1 Tomography principles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 5.4 Coding of a fast inversion routine . . . . . . . . . . . . . . . . . . . . . . . . . . 13

5.4.1 Minimizing the Fisher information . . . . . . . . . . . . . . . . . . . . . 13 5.4.2 Establishment of the good λ . . . . . . . . . . . . . . . . . . . . . . . . . 14

5.5 Postreconstruction test for bad channels . . . . . . . . . . . . . . . . . . . . . . 14

6 Illustrations for shot 21515 15 7 Conclusion 17

7.1 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

A Included files and its history 19

A.1 Necessary files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 A.1.1 fabcat.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 A.1.2 fabsig.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 A.1.3 fabrec.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 A.1.4 Unchanged routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

A.2 Optional files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

A.2.1 Unchanged routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 A.3 Additional tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

B The program code 21

B.1 Necessary files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 B.1.1 fabcat.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 B.1.2 fabsig.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 B.1.3 fabrec.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 B.1.4 fabstore.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

B.2 Optional files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 B.2.1 fabpower.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 B.2.2 fablcfs.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 B.2.3 fabsetetendue.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 B.2.4 fabsetchord.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 B.2.5 fabsetsolidangle.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 B.2.6 fab3d.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 B.2.7 fabproject.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 B.2.8 fabsett.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 B.2.9 fabangle.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 B.2.10 fabstandard.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

B.3 Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 B.3.1 fabmov.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 B.3.2 fabreader.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 B.3.3 fabanalyze.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 B.3.4 fabplot.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

List of Figures

1.1 btomo.m package. Graphical user interface for bolometry tomography reconstruction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2.1 Bolometer design. Legend: 1 - gold-meander resistor; 2 - gold-absorber layer; 3 - golden heat-resistor layer; 4 - kapton carrier foil; 5 - viewing window; 6 - aluminium cooling plate; 7 - golden thermal-contact layer; 8 - conduction path; 9 - mounting support. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.2 Setup of the foil bolometers. Poloidal section of the TCV tokamak. 5 pinhole cameras and the viewing lines of the 64 gold foil bolometers. . . . . . . . . . . . 2

3.1 Main program fabcat.m and its subroutines. The green blocks represent the indispensable subroutines. Yellow units are executed separately by the user to create the variable files (drawn in white) loaded by the main code. . . . . . . 3

3.2 T-matrix calculator. The T-matrix elements T_{ℓi} give the weight of pixel i in

the viewing line of bolometer channel ℓ. . . . . . . . . . . . . . . . . . . . . . . . 4 3.3 Raumwinkel file calculator. The solid angles of each bolometer is calculated. 4 3.4 Bolometric movie player. Permits to check the results for each individual

timeframe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

5.1 Iterations to find the emissivity g. For the minimization of the Fisher information n = 3 loops (flag loopfisher) are enough, for the smoothing loop 4 iterations (flag iterfisher) are executed by default. . . . . . . . . . . . . . . . . 13

6.1 shot 21515 @ t = 0.35 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 6.2 shot 21515 @ t = 0.684 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 6.3 shot 21515 @ t = 0.85 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 6.4 shot 21515 @ t = 1.017 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 6.5 shot 21515 @ t = 1.197 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 6.6 Magnetic flux surface @ 1.197 s. Plasma shot 21515. . . . . . . . . . . . . . 16 6.7 Evolution of the radiated power. Plasma shot 21515. The three curves

represent the total radiated power, the power emitted from the x-point and from the divertor region. Values in [W]. . . . . . . . . . . . . . . . . . . . . . . . . . . 16

6.8 Comparison of the total radiated power and the ohmic heating power. Plasma shot 21515. Values in [W]. . . . . . . . . . . . . . . . . . . . . . . . . . . 16

A.1 Necessary package files for fabcat.m. Located in the main directory \ fabcat. 19 A.2 Files to recalculate the solid angle files of the bolometers. Located in

the subdirectory \ fabcat \ fabsolidangle. . . . . . . . . . . . . . . . . . . . . . . 20 A.3 Fast bolometry result movie player. Located in the subdirectory

\fabcat\fabmovie. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 LIST OF FIGURES iv

A.4 Additional Tools. Located in the subdirectory \fabcat\fabtools. . . . . . . . . 20

Chapter 1

Introduction

The Fast-Algoritm Bolometric Computer Aided Tomography (FABCAT^{1}) reconstruction code is a modified release of the btomo.m package [1] without GUI interface and using a fast minimum Fisher regularization algorithm. The btomo.m package (see figure 1.1) is a useful tool for tomography reconstruction, but it needs quite a lot of computing time and user interaction before delivering results.

Figure 1.1: btomo.m package. Graphical user interface for bolometry tomography reconstruction.

In order to evaluate a plasma shot shortly after the acquisition, one needs a faster tool for the plasma shape reconstruction.

Recently a fast tomography reconstruction algorithm for Soft X-Ray tomography has been developed (the cattcv.m package [2]), but an adequate release for bolometry was still missing.

The new fabcat.m package is simply calling an adaption to the requirements of bolometry of the fast Minimum Fisher Regularization routine from cattcv.m. It also uses an averaged, single reconstruction matrix M which can be applied for all timeframes, making the inversion by the number of timeframes faster (about three seconds in the code at issue).

After testing the package for some particular shots we processed shots 19000 up to 21739 and stored their reconstructions in the MDS tree. In this interval there are 1589 shots, and 668 (42 %) have been reconstructed successfully.

The next chapter gives a short overview and memento on bolometry. Chapter 3 describes the parts and organization of the package for a first orientation. For most of the users of the package the chapter 4 giving operating instructions will probably be the most important. Chapter 5 gives more details on the code, mainly it explains the tomography and their fast formulation for bolometry. Chapter 6 shows reconstructed images for selected shots. For those already familiar with other tomography packages (btomo.m for bolometry, cattcv.m for Soft X-Ray) the chapter A explaining the roots of their codes and its modification will be relevant.

Finally the appendix lists the program code.

This documentation is available online at http://crpplocal/diags/fabcat.pdf (Adobe Acrobat file). The package is administrated by Jan Mlynář and ready to use on HAL under \home\mlynar\matlab\fabcat\.

1colloquial: ’fabulous cat’

Bolometry on Tokamak à configuration variable (TCV)

Bolometric measurements on tokamaks are useful for getting the distribution of the sources of radiation inside the vessel. The TCV at CRPP has a major radius of R = 88 cm and an almost rectangular cross section with width and height of r = 55 cm and z = 155 cm respectively. 8 bolometer arrays with 8 channels each (see figure 2.1) are grouped around a poloidal section of the tokamak and arranged in arrays of 5 pinhole cameras as shown in figure 2.2. Each camera measures the line integrated radiation. A tomographic inversion of the radiation signals provides a spatial resolution of the radiation sources (the emissivity profile).

Figure 2.1: Bolometer design. Legend: 1 - goldmeander resistor; 2 - gold-absorber layer; 3 - golden heat-resistor layer; 4 - kapton carrier foil; 5 - viewing window; 6 - aluminium cooling plate; 7 - golden thermal-contact layer; 8 - conduction path; 9 - mounting support.

The line integrated energy flux is derived from the temperature rise and thermal capacity of well-defined golden foils. However, their thermal inertia limits the temporal resolution of this type of bolometer to approximately 10 ms. The measured signal corresponds not only to the radiated energy, but also contains the contribution due to neutral particle flow [[3]].

The acquisition electronics works at a frequency of 50 kHz and the signals are stored in the MDS tree.

Figure 2.2: Setup of the foil bolometers. Poloidal section of the TCV tokamak. 5 pinhole cameras and the viewing lines of the 64 gold foil bolometers.

Organization of the code

The algorithm is based on the btomo package. Most of the parts has been recycled and adapted, mainly by eliminating the GUI related parts of the code. The fast inversion

Figure 3.1: Main program fabcat.m and its subroutines. The green blocks represent the indispensable subroutines. Yellow units are executed separately by the user to create the variable files (drawn in white) loaded by the main code.

code based on minimizing the Fisher information is essentially the same as that used in Soft-X Ray tomography [2].

3.1 Main program

Figure 3.1 gives an overview of the structure of the code by and names its subroutines.

fabcat.m is the main program and should be called up as described in chapter 4.1 (see paragraph 5.1 for further description).

The subroutine fabsig.m gets the measured bolometry signals from the MDS database, checks for bad channels and does data preprocessing, i.e. corrections respecting the thermal inertia of the bolometers (see paragraph 5.2 for further description).

fabrec.m calculates the tomography inversion using a time averaged M-matrix based on minimization of the Fisher information as regularization functional (see paragraph 5.3 for further description).

fabstore.m stores the reconstruction image (the emissivity g) and the total and partial (from the x-point and divertor) power radiation in the MDS tree. (see paragraph 4.3 about reading the MDS tree).

The .mat files fabchord (coordinates of the viewing lines of the bolometry detectors), fabt (contains the T-matrix) and fabetendue (contains the etendue of the detectors for bolometry) are loaded by the main program and can be recalculated using the associated .m files (see the following paragraph). The precalculated T-matrix was established for a rectangular grid of n_{pixel} = 280 pixel^{1}.

3.2 Supporting routines

3.2.1 fabsetchord.m

Calculates the coordinates of the viewing lines of the bolometers and stores them to fabchord.mat.

3.2.2 fabsetetendue.m

Calculates the bolometers etendue and stores them to fabetendue.mat. The etendue is used to calculate the errorbars of the signals.

3.2.3 fabsett.m

Generates the T-matrix for the inversion, based on the solid angles and etendue of the bolometers. The program loads precalculated fabsolidangle #.mat files and creates the fabt.mat file. See figure 3.2 for the programs hierarchy.

Figure 3.2: T-matrix calculator. The T-matrix elements T_{ℓi} give the weight of pixel i in the viewing line of bolometer channel ℓ.

The program calls the fabsolidangle.m file, which is exactly the same as in the btomo.m package (btt omgrid.m). The fabsolidangle #.mat files can be recalculated using fabsetsolidangle.m as presented in figure 3.3.

Figure 3.3: Raumwinkel file calculator. The solid angles of each bolometer is calculated.

3.2.4 Bolometry movie player

The movie player fabmov.m shows a movie for a given shot and gives the values of χ^{2}, λ and radiated power for each timeframe. It also notes which channels have been deleted.

Figure 3.4: Bolometric movie player. Permits to check the results for each individual timeframe.

1n_{x} = 10 in horizontal and n_{y} = 28 in vertical direction.

3.2.5 Other tools

fabreader.m

The program fabreader.m extracts the status for a shot or for an interval of shots from the MDS tree. Information about the existence of the shot, its successful reconstruction and if the x-point determination was possible is stored in a file fabresult.mat.

fabanalyze.m

The program fabanalyze.m can be used to analyze the results provided by fabreader.m. For the shot interval 19000 to 2173 there are 1589 shots, for 816 of them a x-point was determined and 668 have been reconstructed successfully.

fabplot.m

This tool plots time evolution of plasma current, plasma density and radiated power (total power, power from the X-point and divertor region) read from MDS (prior calculation by fabcat.m is necessary).

Operating instructions

4.1 How to execute ’fabcat.m’

The program can be run in three modes:

- fabcat(shot):

treats the specified plasma shot and delivers the reconstruction at the LIUQE times. The reconstructed image (dimension n_{x} × n_{y} × timeframes, unity is W/m^{3}) is available through the vector X_{inv} (matlab variable name).

- fabcat(shot, t_{min}, t_{max}):

loads and treats each timeframe located between t_{min} and t_{max} [s] only.

- fabcat(shot, t_{min}, t_{max}, points): loads and treats the timeframes located between t_{min} and t_{max} [s] and delivers the results in points frames only. Of course points has to be smaller than the number of acquired images between t_{min} and t_{max}.

4.2 How to configure ’fabcat.m’

In the header of the main program file fabcat.m several flags are specified to control the behavior of the package. When experiencing problems with the execution of the package the first thing to check are these configuration flags and try to run them using the default settings. When you want to change the default settings copy all files to a local directory and run it from there.

- f cam Specifies which pinhole cameras should be used. Each camera gathers eight bolometers looking through the same pinhole of the section (see chapter 2). Default setting = [1 1 1 1 1 1 1 1], i.e. all eight cameras are activated.

- chantest The package can itself try to identify bad channels. The following values are admitted : 0: no automatic testing for bad channels. 1: Testing using Christian Deschenaux’s method (backfiltering). 2: Testing using Arno Refke’s method (consistency check). 3: Testing using both methods. Default setting = 2. For a description of the testing methods see paragraph 5.2.5. As a new feature, the subroutine fabrec.m is also testing for faulty channels (see paragraph 5.5, the configuration is done using the flag LIMIT).

- smoothmode Specifies the method by which the raw bolometry signals are treated and calibrated. 0: no data treatment is done. 1: Processing using Bernard Joye’s method (polynomial fitting). 2: Processing using Christian Deschenaux’s method (backfiltering).

Default setting = 2. For a detailed description of the data processing methods see paragraph 5.2.6.

- delchanuser Allows the specification of channel numbers which will be ignored by the reconstruction. Maybe useful if single bolometers are defective. Default setting = [], i.e. no channel excluded.

- runsigma Provides the global multiplication factor for the error bars. The constants unity is in percent. The variance σ_{i} will still be divided by the individual channel etendue (see formula 5.1). Default setting = 3.

- lambda init Provides the initial λ for the fast minimum Fisher regularization routine, i.e. the balance between fitting and smoothing (see paragraph 5.4.2). Default setting = 0.1.

- lamlim Minimum value of the final λ which qualifies the shot still as a good one. This parameter was chosen arbitrary. Default setting = 0.01.

- chitar Target χ^{2}, i.e. the reconstruction tries to obtain this final χ_{target}^{2}. Default setting = 1, i.e. a reconstruction image having maximum variance σ_{i}.

- errchi Supreme difference between the time averaged final χ^{2} and the target χ_{target}^{2} (see flag chitar) which qualifies the shot still as a good one. This parameter was chosen arbitrary. Default setting = 0.05.

- loopfisher Provides the number of times the main loop of fabrec.m minimizing the Fisher information is going to be executed. Default setting = 3.

- iterfisher Provides the number of times the subloop of fabrec.m (finding λ which brings χ^{2} next to χ_{target}^{2}) is going to be executed. Default setting = 4.

- limit Limiting factor for the postreconstruction bad channel check in fabrec.m. See paragraph 5.5 for its definition. Default setting = 5.

- i bord If set to 1 then emissivity from borders of pixels mesh are constrained to be zero. Although this feature work well for Soft X-ray, it is not suitable in the case of bolometry. Default setting = 0 (disabled).

- x blind If set to 1 then the radiation originating from the X-point region is subtracted from the bolometry signals and a second reconstruction is initiated. Please note that this feature depends on geometry. This part has to be changed if the pixel setup is modified. See the associated subflag cutoff for further configuration. Default setting = 0 (disabled).

- cutoff Specifies the expanded X-point geometry: 0: the X-point is found in the left center and diverter region and pixels radiating down to the percentage defined by the flag level of the maximum radiation are assigned to the expanded X-point. 1: pixel 12 and 13 as expanded X-point. 2: pixels 10 to 16 as expanded X-point. 3: pixels 6 to 18 and 41 as expanded Xpoint. 4: pixles 5 to 20 and 40 to 43 as expanded X-point. Default setting = 0 (automatic determination).

- level Pixels radiating down to the percentage given by this flag of the maximum radiation are assigned to the expanded Xpoint. Default setting = 80 [%].

- radiation Calls the slow fabpower.m (btboloti pradinout.m in the btomo.m package) to calculate the radiated power in regard to the LCFS. Default setting = 0 (disabled).

- saving If set to 1 then stores shotnumber, the emissivity g (Xinv) and the timeframe vector (tvec) to a .mat file on the hard disk. Please specify the path where the files should be stored in using the flag savepath. Default setting = 0 (disabled).

- MDS If set to 1 then the results are stored in the MDS tree. Make sure that you have MDS writing privileges. Default setting = 1 (enabled).

- plotting If set to 1 then makes a contour plot for the last timeframe of the reconstructed emissivity distribution on the screen. Default setting = 0 (disabled).

4.3 How to read the MDS tree

The subroutine fabstore.m (taken from the btomo.m package) stores the data in the MDS tree. The results can be read from MDS using the following nodes:

- \\ results :: btomo : method Information about the used method. There are three possibilities:

- ”FABCAT with X-point at pixel 11”. The X-point has been found at pixel 11, but its radiation emissivity has not been excluded during the inversion.

- ”FABCAT without X-point at pixel 11”. The X-point has been found at pixel 11 and its radiation emissivity has been excluded during the inversion.

- ”FABCAT, X-point determination failed”. No X-point has been determined.

- \\ results :: btomo : comment Gives the final λ and the mean value of χ^{2}. If the shot was reconstructed successfully, a statement e.g. ”shot 15211 OK” is included.

- \\ results :: btomo : rmesh. Pixels center coordinates in radial direction. Dimension n_{x}. Values in [m].

- \\ results :: btomo : zmesh. Pixels center coordinates in vertical direction. Dimension n_{y}. Values in [m].

- \\ results :: btomo : lambda Final smoothing factor λ, see paragraph 5.4.2 for comments.

- \\ results :: btomo : sigma The global multiplication factor for the error bars, as defined by the flag runsigma in the main program (see chapter 4.1). Value in percent.

- \\ results :: btomo : chi squared The final fitting coefficients χ^{2}. If the individual (i.e. for each timeframe) χ^{2} are not close to zero, you may rerun fabcat.m for a restricted time interval or verify with the btomo.m package. See paragraph 5.3 for more comments. Dimension timeframes.

- \\ results :: btomo : nchord Vector listing the enabled channels.

- \\ results :: btomo : smoothing Number of rescaled timeframes, only available when the variable points was specified within fabcat.m call.

- \\ results :: btomo : time Timeframe times in [s].

- \\ results :: btomo : emissivity The emissivity matrix g (Xinv). Dimension n_{pixel} × timeframes. Values in [W/m^{3}]. Notice that this node can be shortened to \\ results :: btomo for future shots to conform with ANASRV standards.

- \\ results :: btomo : prad above At present radiation power coming from the X-point (only saved if the determination of the X-point was successful). The power of the pixel with highest emissivity and all neighboring pixels is summed up. Values in [W]. Dimension timeframes.

- \\ results :: btomo : prad below At present radiation power coming from the divertor region below the X-point (only saved if the determination of the X-point was successful). The two leftmost pixel columns are summed up from the bottom of the vessel up to two pixels below the pixel of the x-point with highest emissivity. Values in [W]. Dimension timeframes.

- \\ results :: btomo : prad tot Total radiation power of the whole tokamak. Values in [W]. Dimension timeframes.

Detailed description of the method

5.1 fabcat.m

The header of the main program list all the necessary user adjustable configuration flags. Their signification has been described in chapter 4.

The MDS data is read by fabsig.m (see next paragraph) for the specified shot with respect to the chosen cameras, channels and times.

The precalculated T-matrix (including the pixel grid, all values in [cm]) is stored in fabt.mat and simply loaded. The error bars df for the inversion are calculated using the etendue (A · Ω) of the bolometry detectors (stored in fabetendue.mat):

df = σ · \frac{sup(A · Ω)}{(A · Ω)}  (5.1)

The inversion is done using the subroutine fabrec.m (see paragraph 5.3).

The resulting tomography reconstruction (X_{inv}) is converted to [W/m^{3}]. χ^{2} and λ are tested for the quality^{1} of the reconstruction (by default^{2} a good quality image has λ > 0.01 and | < χ^{2} > − 1 | < 0.05) and commented on screen.

According to the configuration flag plotting a contour plot of X_{inv} (see figure 6.1 ff.) can be displayed or not.

The power is calculated from the emissivity by the integration over the ring tubes intersected by the pixels. The power coming from a rectangular shaped ring (cylindrical geometry, cross section dr·dz) intersected by pixel positioned in (r, z) is

P (r, z) = ∫_{0}^{2π} ∫_{left}^{right} ∫_{low}^{up} g(r, z)·r·dθ·dr·dz

(5.2)

where θ is the toroidal angle which gives the position on the rings central line. With ’right’ we mean the right border of the pixel. Assuming that the power emission is constant along the rings central line (cylindrical symmetry) and with pixels having horizonal and vertical extension of dx and dy respectively leads to

P (pixel) = 2π · g(pixel) · r(pixel) · dx · dy (5.3)

Summing up all n_{pixel} pixels provides the total radiated power P_{tot}:

P_{tot} = ∑_{pixel=1}^{n_{pixel}} P (pixel)  (5.4)

5.2 fabsig.m

This subroutine is responsible for extracting and preparing any information from the MDS tree. It is the only routine using MDS reading access commands.

1i.e. the mean χ^{2}, averaged over all timeframes.

2see the empirical parameters errchi and lamlim as described in chapter 4

The behavior of fabsig.m is essentially determined by the number of parameters provided when executing fabcat.m (see chapter 4.1).

5.2.1 The brightness f_{ℓ}

During a shot, the plasma radiates. The radiation emitted in r⃗ is characterized by the emissivity^{3} g(r⃗).

Detector ℓ collects the total radiated power P_{ℓ} from its field of view given by the solid angle Ω_{ℓ}. Hence, the total radiated power can be expressed as an integral of the cone of view S_{ℓ} defined by the detector and aperture geometry

P_{ℓ} = ∫_{S_{ℓ}} dr⃗ g(r⃗) · \frac{Ω_{ℓ}(r⃗)}{4π}, ℓ = 1, .., 64 (5.5)

The bolometers have a narrow field of view. The emissivity g(r⃗) doesn’t vary considerable on a surface perpendicular to the axis^{4} of the cone S_{ℓ}. The volume integral can be replaced by a simple line integral dr⃗ → A(s) ds where A is the surface of the detector and ds is a line element along the line of sight s_{ℓ}.

Using the brilliancy theorem^{5} for the etendue (A · Ω) we can write

P_{ℓ} = \frac{(AΩ)_{ℓ}}{4π} ∫_{s_{ℓ}} ds g(r⃗), ℓ = 1, .., 64 (5.6)

The multiplication with factor 1.78 · 10^{4} in the code corresponds to the determination of the chord brightness f_{ℓ} which is derived from the total radiated power P_{ℓ} as follows

f_{ℓ} = \frac{P_{ℓ}}{(AΩ)_{ℓ}/4π} , ℓ = 1, .., 64  (5.7)

In the following the measured chord brightness f will be called signals.

5.2.2 Reading the MDS tree

The following MDS trees are read :

- \\ base :: bolo : signals The matrix (size channels × timeframes) contains preprocessed bolometry signals, i.e. the acquired chord brightness f_{ℓ} integrated along the lines of sight of the detector ℓ.

- \\ base :: bolo : source The raw bolometer signals f̆_{ℓ} matrix (size channels×timeframes). They are used for weird channel detection.

- dim of(\\ base :: bolo : signals) The normal time base vector (length timeframes) containing the the times corresponding to the timeframes of the signals matrix.

- dim of(\\ results :: r axis) The LIUQE time base vector. These are the times used for the magnetic field reconstruction of the LIUQE package.

- \\ base :: bolo : tau τ, the time constants of the bolometers. They are used for data treatment.

5.2.3 LIUQE times

When fabcat.m has been called providing the shot number only, then the LIUQE times are used for the reconstruction. A new time base consisting of normal times nearest to LIUQE times is calculated.

5.2.4 NORMAL times

When fabcat.m is called providing a minimum time argument (t_{min}) then the normal times are used. If no maximum time (t_{max}) was specified all timeframes beginning with t_{min} are loaded. The optional argument points is to use only points timeframes for the specified shot. This decreases calculation time. Its clear that points has to be lower than the total timeframes available for a given shot.

3power per volume

4the so-called line of sight

5A′ · Ω′ = A · Ω through an optical system

5.2.5 Automatic channel deleting

The header of the main program fabcat.m allows for the specification of channels which should not be used for reconstruction. Using the flag channels in the main program (see chapter 4.1) it is possible to choose all, none or one of the following automatic channel analysis algorithm to identify further faulty channels:

- Arno Refke’s method. The raw raw signals are checked for consistency, i.e. the signal range, signal mean value, signal slope and negative signals.

- Christian Deschenaux’s method. Filters the raw signals using a butterworth 4^{th} order digital lowpass filter and checks the backfiltered signals mean values for consistency.

These channel deletion methods are the same as those used in the btomo.m package.

5.2.6 Signal processing

The flag smoothmode (in the header of fabcat.m) can be used to choice if one, none or one of the following smoothing procedures is used to calibrate the measured signals:

- Bernard Joye’s method. Fits the temporal behavior of each channel with a polynomial of 2^{nd} degree over 30 time samples considering the bolometers time constant τ for the linear term.

f_{ℓ} = at^{2} + b\frac{τ}{△t}t

where a is the quadratic and b the linear fit coefficient. △t is the averaged time step per frame.

- Christian Deschenaux’s method. Filters the raw signals again using a butterworth 4^{th} order digital lowpass filter :

f_{ℓ} = f̆_{ℓ} + τ · \frac{df̆_{ℓ}}{dt}

where f̆_{ℓ} designates the backfiltered raw signals.

5.3 catrec.m

This is the core piece of software of this package. The inversion is done under the control of several parameters specified in the main program and described in chapter 4.1. The routine corresponds to the one used in the Soft X-ray package and has only been modified slightly (see appendix A).

5.3.1 Tomography principles

The tomography problem consist of solving the system of integral equations

f_{ℓ} = ∫_{s_{ℓ}} ds g(r⃗), ℓ = 1, .., 64  (5.8)

where f_{ℓ} are known and g(r⃗) is unknown.

The system 5.8 of inhomogeneous Fredholm equations of the first kind [3] is always underdetermined, since an infinite number of line integrated data measurements f_{ℓ} would be necessary to determine g(r⃗) exactly. By a subdivision of the area of the poloidal cross section into pixels, the number of degree of freedom is reduced. We use a grid of rectangular pixels (by default n_{x} = 10 pixels in radial and n_{y} = 28 in vertical direction) covering the whole vessel of the tokamak. The system 5.8 has now been transformed to a system of algebraic equations^{6}

f_{ℓ} = T_{ℓi} · g_{i}, ℓ = 1, .., 64  (5.9)

where T_{ℓi} describes the weight of pixel i in the line of sight of detector ℓ. The Tmatrix is related to geometry only (detector positions, aperture dimensions and pixel grid) and is calculated using fabsett.m (see paragraph 3.2.3).

g represents now the two-dimensional emissivity distribution.

6sum over repeated indexes (Einsteins convention).

Unfortunately a direct inversion of T is not possible because it is badly conditioned. Furthermore, with our grid we have n_{pixel} = 280 > n_{ℓ} = 64, so the set of equations 5.9 is overdetermined and none or an infinite number of solutions is possible. To get a unique and sensible solution of the problem a functional

ϕ = \frac{1}{2}χ^{2} + λℜ  (5.10)

is minimized (see [4]). χ^{2} is calculated by

χ^{2} = ∑_{i} [(\frac{T_{ℓi} · g_{i} − f_{i}}{σ_{i}} )^{T} · (\frac{T_{ℓi} · g_{i} − f_{i}}{σ_{i}} )]

(5.11)

and controls the goodness^{7} of fit. ℜ is a regularization functional controlling the smoothness of the fit. The smoothing parameter λ weights between goodness of the fit and smoothness and has to be chosen carefully for reliable results.

The ℜ functional in fabrec.m is the minimum of the Fisher Information I_{F}

I_{F} = ∫ dx\frac{(g′(x))^{2}}{g(x)}  (5.12)

leading to the smoothest [5] solution of the inversion. The division by g(x) in 5.12 assures that regions of low g are particularly smooth, taking care of low-emissivity regions.

5.4 Coding of a fast inversion routine

The process to establish the emissivity g was applied to each timeframe t in the btomo.m package. The essential reduction in time consumption of fabrec.m is related to averaging χ^{2} and λ and the reconstruction matrix M over the whole shot. The routine uses two nested loops which are shown in figure 5.1.

Figure 5.1: Iterations to find the emissivity g. For the minimization of the Fisher information n = 3 loops (flag loopfisher) are enough, for the smoothing loop 4 iterations (flag iterfisher) are executed by default.

5.4.1 Minimizing the Fisher information

The algorithm uses a variational principle to minimize the Fisher information. (g′(x))^{2} = (\frac{d}{dx}g(x))^{2} in 5.12 can be evaluated as gradients in respect to x and y. Assuming that ∇_{x} and ∇_{y} are finite-difference matrix representations of the partial derivatives^{8}, then

(g′(x))^{2} = (∇_{x}g)^{T} ∗ (∇_{x}g) +
(∇_{y}g)^{T} ∗ (∇_{y}g) (5.13)

The weight factor \frac{1}{g(x)} in 5.12 can be introduced using a diagonal weight matrix W.

But to avoid nonlinearity we start with

W_{ij}^{(1)} = δ_{ij}  (5.14)

and use the found g^{(1)} (found after the first iteration) and so on, i.e.

W_{ij}^{(n)} = \frac{1}{G_{i}^{(n−1)}} · δ_{ij}, n ⩾ 2  (5.15)

where

G^{(n)} = ∑_{t} g^{(n)}(t)  (5.16)

is the over all sum of the timeframes of vectors of emissivity g in a shot. By means of

H^{(n)} = (∇_{x})^{T} ∗ W^{(n)} ∗ (∇_{x}) +

(∇_{y})^{T} ∗ W^{(n)} ∗ (∇_{y}), n ∈ ℕ (5.17)

7low χ^{2} corresponds to ’overfitting’, χ^{2} = 1 is the solution with maximum variance.

8Bx and By in the code

5.10 can now be rewritten as

ϕ^{(n)} =

\frac{1}{2} (\frac{T∗G^{(n+1)}−F}{σ} )^{T} ∗ (\frac{T∗G^{(n+1)}−F}{σ} ) +

λ (G^{(n+1)}^{T} ∗ H^{(n)} ∗ G^{(n+1)})  (5.18)

where

F = ∑_{t} f(t)  (5.19)

is the over all sum of the timeframes of vectors of the chord brightness f. Minimizing ϕ^{(n)} means

\frac{∂ϕ^{(n)}}{∂G_{i}} = 0 ∀ i = 1, .., 280  (5.20)

leading to the normal equations

(\frac{T^{T} ∗ T}{σ^{2}} + λH^{(n)}) ∗ G^{(n+1)} = \frac{T^{T} ∗ F}{σ^{2}}

(5.21)

which has to be solved iteratively for G^{(n+1)}. In fabrec.m [2] this is done by calculating a reconstruction matrix M (\ stands for left division)

M^{(n)} = (\frac{T^{T} ∗ T}{σ^{2}} + λH^{(n)}) \ \frac{T^{T}}{σ} (5.22)

The shot-averaged reconstruction matrix M can now simply be applied to the single timeframe

g^{(n)}(t) = M^{(n)} ∗ f(t)  (5.23)

5.4.2 Establishment of the good λ

Before the next loop minimizing I_{F} is executed (using a new weight matrix W determined by g^{(n−1)}) a new λ which iterates χ^{2} to χ_{target}^{2} has to be found.

In general a solution with maximum variance σ requires χ^{2} = 1. But due to our definition of χ^{2} (see formula 5.11) we have summed over the n_{ℓ} enabled bolometer channels (by multiplying the matrices) and also over the n_{t} timeframes of the shot (by time averaging)

χ_{target}^{2} = n_{ℓ} · n_{t}  (5.24)

For a given iteration n of the invoking Fisher loop the m^{th} subloop simply calculates

λ^{(m+1)} = λ^{(m)} + \frac{dλ^{(m)}}{dχ^{2}} (χ_{target}^{2} − χ^{2}) (5.25)

with

\frac{dλ}{dχ^{2}} = \frac{λ^{(m−1)} − λ^{(m−2)}}{χ^{2}^{(m−1)} − χ^{2}^{(m−2)}}  (5.26)

and

λ^{(n=1,m=1)} = λ_{initial} · \frac{T r(T)}{T r(H)}

λ^{(n,m=2)} = \frac{1}{2}λ^{(n,m=1)} (5.27)

Of course evaluating χ^{2} requires the recalculation of the new g^{(n,m)}(t) with the new λ^{(m)} as described in the previous paragraph.

5.5 Postreconstruction test for bad channels

Due to problems with the weird channel detection algorithm described in chapter 5.2.5, we modified the fabrec.m routine in order to test for weird channels after the reconstruction, i.e. by comparing the measured signals and reconstructed emissivity. When we define this difference by

χ = T ∗ g − f  (5.28)

then a channel is declared as bad using the criterion

|m(χ)| ⩾ limit · s(m(χ))  (5.29)

where m(χ) is the mean value over all timeframes and s(m(χ)) its standard deviation for a given channel. By default the parameter limit is set to 5.

Illustrations for shot 21515

In this chapter some results for shot 21515 are presented.

The left border of the cross section is oriented to the center of the TCV tokamak.

The figures 6.2 up to 6.5 shows the evolution of the emissivity g for the shot 21515. The colorbar of these figures are all scaled identically. The values are given in [W/m^{3}]. The images have been reconstructed using the default settings of the package.

Figure 6.1: shot 21515 @ t = 0.35 s.

Figure 6.2: shot 21515 @ t = 0.684 s.

Figure 6.3: shot 21515 @ t = 0.85 s.

Figure 6.4: shot 21515 @ t = 1.017 s.

Figure 6.5: shot 21515 @ t = 1.197 s.

CHAPTER 6. ILLUSTRATIONS FOR SHOT 21515 16

Figure 6.6 shows the magnetic configuration of the plasma shot at t = 1.197 s. Figure 6.7 shows the evolution of the radiated power for the whole vessel, the x-point and from the divertor region. Figure 6.8 compares the resistive heating power with the total radiated power.

Figure 6.6: Magnetic flux surface @ 1.197 s.

Plasma shot 21515.

Figure 6.7: Evolution of the radiated power. Plasma shot 21515. The three curves represent the total radiated power, the power emitted from the x-point and from the divertor region. Values in [W].

Figure 6.8: Comparison of the total radiated power and the ohmic heating power. Plasma shot 21515. Values in [W].

Chapter 7 Conclusion

The new package fabcat.m has been realized and tested for a large number of shots.

Suitable default parameters were chosen. The results can be stored to the MDS database and simply visualized, so that the package is ready for daily use at CRPP.

The structure of the package is simple enough and its routines have been well commented so that it is ready for future development.

To make the package still more useful, it should be integrated on ANASRV for automatic post-shot processing. Traces of radiated power could then make part of TCV operators scope.

Studies on the latest shots have shown that fabcat.m succeeds in about 40 % of the shots to deliver a good quality image reconstruction. Random spot checks have shown that even the btomo.m package fails at the same shots. This is due to bad signal data, i.e. systematic problems with the top and bottom bolometer arrays. Regular calibration of the bolometers is necessary, but those above and below the vessel are not easily accessible.

In future we would like to include the fluxshaped pixels method described in [7]. But detailed studies of this novel method within btomo.m (where this is already implemented) is necessary first.

7.1 Acknowledgements

We would like to thank the TCV group, namely Richard Pitts, Bernard Joye and Basil Duval for their support. Thanks also to all them who contributed with programming to previous packages solving the tomography problem in plasma physics.

Christian Schlatter and Jan Mlynář

Lausanne, february 2002

References

[1] Iwan Jerjen and Jan Mlynář; BTOMO: The Software Package Used for the Bolometry on TCV Tokamak, CRPP EPFL Lausanne, INT 197/99.

[2] Jan Mlynář; Rapid tomography inversion and its contribution to plasma position measurements at the TCV tokamak; not published yet.

[3] Jan Mlynář, Ivo Furno, Bernard Joye and Arno Refke: Bolometry on the TCV Tokamak, CRPP EPFL Lausanne, INT 196/99.

[4] Martin Anton and al.:X-ray tomography on the TCV tokamak, Plasma Physics Controlled Fusion, Volume 38, 1996.

[5] B. R. Frieden: Application to optics and wave mechanics of the criterion of maximum Cramer-Rao bound, Journal of Modern Optics, Volume 35, page 1297, 1988.

[6] William H. Press and al:Numerical recipes in FORTRAN 90, Cambridge 1996.

[7] S Kalvin and S Zoletnik: Tomographic Reconstruction of Radiation Intensity Profiles from Bolometric Measurements on the TCV Tokamak, CAT-Science Bt. for CRPP EPFL Lausanne, INT.

Appendix A

Included files and its history

This chapter list all the files included in the fabcat.m package. It is also directed at those interested in the history of the code. Here we will explain the modifications of already existing routines (packages btomo.m and cattcv.m) to meet the requirements of a fast bolometry inversion.

A.1 Necessary files

Figure A.1: Necessary package files for fabcat.m. Located in the main directory \ fabcat.

Figure A.1 shows the necessary files to run the package. If the .mat files (the files without an extension in figure A.1) don’t exist already, they can be created using the associated .m files.

A.1.1 fabcat.m

The main program is fundamentally based on the main program of the cattcv.m package. All lines of code related to Soft-X ray tomography only have been eliminated. Some elementary setup parameters located in catsetup.m for the cattcv.m package have been moved to the main program. So an additional setup routine is no longer needed. New, added features are the saving in the MDS tree and the fast calculation of the radiated power.

A.1.2 fabsig.m

Essential reduction of the routine btomodata.m of the btomo.m package. All GUI related elements have been removed. Some minor adaption concerning the data types of the variables were necessary.

A.1.3 fabrec.m

Ancient fastfish.m of the cattcv.m package. Addition of the i bord flag. The zero radiation border can now be disabled. Several modifications were necessary (see all ’if i bord ...’ expressions).

An post-reconstruction χ^{2} test disabling bad channels has been added, since the Deschenaux backfiltering method does not always work satisfactory.

The definition of the border pixels has been corrected.

The plasma position measurement has been removed.

A.1.4 Unchanged routines

The following routines haven’t been modified but renamed:

- btboloti store mds.m ⟶ fabstore.m (to store results to MDS),

- btboloti pradinout.m ⟶ fabpower.m (to calculate the power radiated inside and outside the LCFS) and

- bttcv getlcfs.m ⟶ fablcfs.m (subroutine used by fabpower.m to determine the LCFS).

A.2 Optional files

The files fabsetetendue.m and fabsett.m are used to create the .mat files fabetendue and fabt and are located in the main directory (see figure A.1). fabsett.m needs the subroutines fabangle.m, fabstandard.m and the 64 files fabsolidangle #.mat from the subdirectory \fabsolidangle. The latter can be created using the fabsetsolidangle.m program in this subdirectory (see figure A.2).

The files fab3d.m, fabproject.m and fabchord are necessary to run fabsetsolidangle.m. The .mat file fabchord can be created using fabsetchord.m if it doesn’t exist already.

A.2.1 Unchanged routines

The following routines haven’t been modified but renamed:

- omgrid 3d.m ⟶ fab3d.m and

- projbl.m ⟶ fabproject.m

Figure A.2: Files to recalculate the solid angle files of the bolometers. Located in the subdirectory \ fabcat \ fabsolidangle.

- btt omgrid new.m ⟶

fabsetsolidangle.m,

- omgrid main new.m ⟶

fabsolidangle.m

A.3 Additional tools

In the subdirectories \ fabtools (see figure A.3) and \ fabmovie (see figure A.3) you will find the tools described in chapter 3.2.5.

The fabmov.m program is a modified version of the catmov.m program for soft X-ray.

Figure A.3: Fast bolometry result movie player. Located in the subdirectory \fabcat\fabmovie.

Figure A.4: Additional Tools. Located in the subdirectory \fabcat\fabtools.

Appendix B The program code

B.1 Necessary files

B.1.1 fabcat.m

function [ Xinv , tvec , xmesh , ymesh , chi2 , lambda , …
FastPtot , FastPbelow , FastPxpoint , good ] = …
fabcat ( shot , tmin , tmax , points )
% f u n c t i o n [ Xinv , tvec , xmesh , ymesh , chi2 , lambda , …
%  FastPtot , FastPbelow , FastPxpoint , good ] = …
%  fabcat ( shot , tmin , tmax , points )

%

% Algorithm f o r f a s t tomography r e c o n s t r u c t i o n f o r bolometry .
% Package Fast−Algorithm Bolometric Computer Aided Tomography (FABCAT) .
% Only one i n v e r s i o n matrix M i s c a l c u l a t e d per shot .
% I f the tmin , tmax times were not s p e c i f i e d , a l l LIUQE times are applied .
% Otherwise a l l ( normal ) bolometry time samples within tmin , tmax are
% r e c o n s t r u c t e d . Points l i m i t s the number of timeframes in the r e c o n s t r u c t i o n .

%

% D e f i n i t i o n of the X−point :
% The p i x e l with h i g h e s t e m i s s i v i t y and i t s upper , lower and r i g h t adjacent
% p i x e l s within p i x e l s 3 to 1 6 and 3 1 to 4 4 ( the i n f e r i o r part of two l e f t
% border columns , the so−c a l l e d p o t e n t i a l X−point r e g i o n ) i s are searched
% and excluded .
% ATTENTION : THESE DEFINITIONS DEPEND ON THE CHOSEN PIXEL GRID ! ! !

%

% See f i l e FABCAT. pdf f o r f u r t h e r documentation ( a v a i l a b l e on http :// c r p p l o c a l ) .

%

% The r e s u l t s can be reached via g l o b a l v a r i a b l e s :

%

% Y  : the processed ( c a l i b r a t e d ) bolometry s i g n a l s
% tvec  : the vector of times
% Xinv  : 3−dim e m i s s i v i t y ( ny x nx x tvec ) r e c o n s t r u c t i o n image [W/mˆ 3]
% xmesh , ymesh : g e o m e t r i c a l setup of the p i x e l s [ cm ]
% chi2 , lambda : r e s u l t i n g chi ˆ 2 and smoothing f a c t o r
%  length ( chi2 )= length ( tvec )

% FastPtot  : Radiated power i n s i d e the whole tokamak [W]
% FastPxpoint : Radiated power coming from the x−point [W]
% FastPbelow  : Radiated power coming from the d i v e r t o r r e g i o n [W]
%  All radiated power v a r i a b l e s have length ( tvec )
% good  : Flag i n d i c a t i n g i f the r e c o n s t r u c t i o n was s u c c e s s f u l

%

% History :

%

% C h r i s t i a n S c h l a t t e r , TPIV CRPP EPFL , February 2 0 0 2 ( f i r s t r e l e a s e )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
global Y tvec Xinv chi2 lambda xmesh ymesh e r r f l a g
global FastPtot FastPbelow FastPxpoint good
disp ( ' ' )

disp ( 'Fast-Algorithm Bolometric Computer Aided Tomography (FABCAT).' ) ;
disp ( 'Final release 1.4' )
disp ( 'Latest modification on 2002/02/07.' )

disp ( ' ' ) ;
debut=cputime ;
% User a d j u s t a b l e parameters ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
f cam  =  [ 1 1 1 1 1 1 1 1 ] ;
% Camera switch . 0 = d i s a b l e d ; 1 = enabled ( d e f a u l t ) .
chantest  = 1; % 0 : no t e s t i n g f o r weird channels .
% 1 : t e s t i n g using C h r i s t i a n Deschenaux ’ s method ( d e f a u l t )
%  based on b a c k f i l t e r i n g of the bolometer s i g n a l s .
% 2 : t e s t i n g using Arno Refke ’ s method
%  based on a c o n s i s t e n c y check of the bolometer s i g n a l s .
% 3 : t e s t i n g using both methods .
smoothmode = 2; % 0 : no data p r o c e s s i n g .
% 1 : Bernard Joye ’ s data treatement , based on f i t t i n g of
%  a polynomial of second degree to the s i g n a l s ( slow ) .
% 2 : Data treatment using C h r i s t i a n Deschenaux ’ s method
%  based on b a c k f i l t e r i n g of the raw s i g n a l s ( d e f a u l t ) .
% All data p r o c e s s i n g methods c a l i b r a t e the s i g n a l in r e s p e c t
% to the bolometer thermal i n e r t i a c o e f f i c e n t s TAU.
delchanuser = [ ] ; % User defined channels to exclude from the r e c o n s t r u c t i o n .
runsigma  = 3; % Errorbar ( f o r the s i g n a l s ) m u l t i p l i c a t i o n f a c t o r .
% The f i n a l e r r o r b a r s are s t i l l m u l t i p l i e d by the
% bolometer etendue . Default = 3 %

lambda init = 0 . 1 ; % I n i t i a l lambda f o r the minimum f i s h e r loop ( d e f a u l t = 0 . 1 ) .
lamlim  = . 0 1 ; % Minimum lambda f o r a good shot ( d e f a u l t = 0 . 0 1 ) .
c h i t a r  = 1; % Target value of chi ˆ 2 . Default = 1
e r r c h i  = . 0 5 ; % Supremum of the d i f f e r e n c e between the time averaged f i n a l
% chi ˆ 2 and the t a r g e t chi ˆ 2 defined by the f l a g CHITAR.
% Default = 0.05
l o o p f i s h e r = 3; % Number of loops of the minimum f i s h e r algorithm to
% minimize the Fisher information ( d e f a u l t = 3 ) .
i t e r f i s h e r = 4; % Number of i t e r a t i o n s of the minimum f i s h e r algorithm
% to f i n d the lambda that pushes chi ˆ 2 to c h i t a r
% ( d e f a u l t = 4 ) .
l i m i t  = 5; % l i m i t i n g f a c t o r f o r post−r e c o n s t r u c t i o n channel
% c o n s i s t e n c y check ( d e f a u l t = 5 ) .
i b o r d  = 0; % 0 : Emissivity c o n t r i b u t i o n coming from the v e s s e l s
%  border p i x e l s are taken in c o n s i d e r a t i o n ( d e f a u l t ) .
% 1 : emission form border p i x e l s i s s e t to zero .
x b l i n d  = 0; % 1 : Excludes the expanded X−point from the r e c o n s t r u c t i o n
%  to i n c r e a s e the c o n t r a s t of the r e c o n s t r u c t e d image .
%  For the d e f i n i t i o n of the X−point see f l a g CUTOFF.
% 0 : no X−point e x c l u s i o n ( d e f a u l t ) .
c u t o f f  = 0; % 0 : Automatic mode ( d e f a u l t ) .
%  The X−point i s searched automatically ( see d e s c r i p t i o n
%  in the header of t h i s f i l e ) and f i n a l l y expanded to
%  p i x e l s in the p o t e n t i a l X−point r e g i o n having
%  e m i s s i v i t i e s higher than the percentage defined by
%  the f l a g LEVEL in r e s p e c t to the X−point .
% 1 : Pixel 1 2 and 1 3 as expanded X−point .
% 2 : Pixel 1 0 to 1 6 as expanded X−point .
% 3 : Pixel 6 to 1 8 and 4 1 as expanded X−point .
% 4 : Pixel 5 to 2 0 and 4 0 to 4 3 as expanded X−point .
l e v e l  = 80; % Percentage of the X−points e m i s s i v i t y to e s t a b l i s h the
% expanded X−point (CUTOFF has to be = 0 ) . Default = 8 0 [ % ] .
r a d i a t i o n  = 0; % 1 : The f r a c t i o n s of r a d i a t i o n i n s i d e and o u t s i d e of the
%  LCFS are determined using the slow fabpower .m
%  r o u t i n e ( o r i g i n a l l y the b t b o l o t i p r a d i n o u t .m r o u t i n e ) .
% 0 : using the f a s t algorithm to determine FastPtot ,
%  FastPxpoint and FastPbelow ( see header of t h i s f i l e ) .

saving  = 0; % 1 : Stores the v a r i a b l e s shot , Xinv , tvec , FastPtot ,
%  FastPbelow , FastPxpoint and good to the harddisk in
%  f i l e s c a l l e d FABCAT shotnumber . mat . Use the f l a g
%  SAVEPATH to chose the d e s t i n a t i o n d i r e c t o r y .
% 0 : no saving to harddisk ( d e f a u l t ) .

savepath  =  '/home/tp47/matlab/FBTR/' ;

% Destination f o l d e r f o r the FABCAT shotnumber . mat f i l e s .
% Make sure that you have w r i t i n g p r i v i l i g e s f o r the chosen
% d e s t i n a t i o n f o l d e r . See f l a g SAVING f o r more d e t a i l s .
MDS  = 0; % 1 : The r e s u l t s of the r e c o n s t r u c t i o n are saved to the
%  MDS−t r e e . Make sure you have MDS w r i t i n g p r i v i l i g e s .
% 0 : f e a t u r e d i s a b l e d ( d e f a u l t ) .
p l o t t i n g  = 0; % 1 : shows a contour p l o t of the r e c o n s t r u c t e d e m i s s i v i t y
%  f o r the l a s t timeframe in the shot .
% 0 : no p l o t t i n g ( d e f a u l t ) .
% ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
e r r f l a g = 0;
good  = 0;

try

close ( 'Emissivity matrix for the last timeframe [W/m^3]' ) ;

end

i f nargin == 1  % LIUQE time mode f o r a s i n g l e shot .
tmin=NaN;
tmax=NaN;

end

i f nargin == 0  % LIUQE time mode f o r multiple shots .
tmin=NaN;
tmax=NaN;

i n p u t f i l e=input ( 'Please provide the name of the shot array : ' , 's' ) ;

try

load ( i n p u t f i l e ) ;

disp ( 'The file has successfully been loaded.' ) ;

catch

disp ( 'The file doesn''t exist. Try another filename.' ) ;

return end

i f ˜ exist ( 'shot' ) % Test i f the provided f i l e i s a r e l i a b l e shot array

disp ( 'This is not a shot array.' ) ;
disp ( 'There wasn''t a variable called shot inside.' ) ;

return end

f i r s t s h o t = shot ( 1 ) ;
l a s t s h o t = shot ( length ( shot ) ) ;
clear shot ;

else

f i r s t s h o t = shot ;
l a s t s h o t = shot ;
clear shot ;

end

i f nargin < 3
tmax = NaN;

end

i f nargin < 4  % Normal time mode without r e s c a l i n g of the timeframes .
points=NaN;

end

for shot = f i r s t s h o t : l a s t s h o t
% FIRST STEP : GETTING THE SIGNALS
Y = [ ] ;
tvec = [ ] ;
goodchans = [ ] ;
e r r s t r = [ ] ;
[Y, tvec , e r r f l a g , e r r s t r , goodchans ] = …
f a b s i g ( shot , f cam , delchanuser , smoothmode , …

APPENDIX B. THE PROGRAM CODE 26

tmin , tmax , points , chantest ) ;
% FABSIG.m g e t s the s i g n a l s from the MDS t r e e , checks f o r weird channels and
% does data p r o c e s s i n g ( c a l i b r a t i o n ) .
% I f TMIN, TMAX were not s p e c i f i e d , s e t s TVEC to n e a r e s t LIUQE times .
% Otherwise reads the normal times from the MDS t r e e .
% I f POINTS was submitted , the number or timeframes i s r e s c a l e d to POINTS
% frames .
i f e r r f l a g  % Shows e r r o r message when a problem occured c a l l i n g FABSIG
disp ( '  ' ) ;
disp ( '  ' ) ;
disp ( '  ' ) ;
disp ( e r r s t r ) ;

return end

% SECOND STEP : CHOOSING GEOMETRY
load fabt  % Saved by the o p t i o n a l r o u t i n e FABSETT.m. I n c l u d e s
% NX, NY ( number of p i x e l s ) ; XMIN, XMAX, YMIN, YMAX
% ( extremal p i x e l p o s i t i o n s ) ; DX, DY ( p i x e l e x t e n s i o n s ) ;
% XMESH, YMESH ( p i x e l s c e n t e r c o o r d i n a t e s ) ; NUMDET
% ( camera i n d i c a t o r ) and T ( the T−matrix ) .
[ nt , nl ] = size (Y) ;
% NT: number of timeframes .
% NL : number of enabled channels .
[ usedch , igoodchans , inumdet ] = i n t e r s e c t ( goodchans , numdet ) ;
% Determines the channels which are good and a c t i v a t e d .
% IGOODCHANS and INUMDET give the index
% of the elements which are common .
% THIRD STEP : DETERMING THE ERRORBARS.
runsigma = runsigma / 1 0 0 ;
% change to [ % ] .
dY = runsigma ∗ ones ( 6 4 , 1 ) ;
% C a l c u l a t i o n of DY, the ’ r e l a t i v e ’ e r r o r of the s i g n a l s .
load fabetendue % Contains the etende of the bolometers in srad ∗cmˆ 2.

% Calculated using FABSETETENDUE.m.
angfact = 1 . / etend ’ ;
wdy = angfact . / min( angfact ) ;
% Weighting matrix f o r the e r r o r b a r s .
dY = dY . ∗ wdy ;  % Errorbars .
Y = Y ’ ;  % because of : s i z e ( framedata ) = ( tvec x goodchans )
% but you need s i z e (Y) = ( goodchans x tvec ) .
% FOURTH STEP : EMISSIVITY RECONSTRUCTION.
lambda = [ ] ;
chim = [ ] ;
Xinv = [ ] ;
chi2 = [ ] ;
xpoint = 1 ;  % The X−point hasn ’ ’ t been excluded yet .
[ Xinv , lambda , chim , inumdet ] = …
f a b r e c (Y, dY , xmesh , ymesh , tvec , inumdet , igoodchans , …
T , lambda init , l o o p f i s h e r , i t e r f i s h e r , i b o r d , l i m i t ) ;
% FABREC.m i s doing the i n v e r s i o n . The r e c o n s t r u c t e d image i s d e l i v e r e d
% by Xinv . Final smoothing f a c t o r LAMBDA and chi2 (CHIM) are returned .
% The r o u t i n e i s checking again f o r weird channels using a post−
% r e c o n s t r u c t i o n t e s t on chi ˆ 2 ( see f l a g LIMIT ) .
Xinv=Xinv ∗ 1 0 0 ;  % Conversion to [W/mˆ 3 ] .
i f p l o t t i n g & ˜ x b l i n d
% Shows a contour p l o t of the found e m i s s i v i t y of the
% f i r s t i n v e r s i o n i f no second i n v e r s i o n i s planned .

   figure ( 'Name' , 'Emissivity matrix for the last timeframe [W/m^3]' ) ;

contourf ( xmesh / 1 0 0 , ymesh / 1 0 0 , Xinv ( : , : , length ( tvec ) ) , 5 0 ) , …
shading f l a t , axis equal
colorbar ;

end

% FIFTH STEP : RECONSTRUCTION RELIABILITY TESTING
chi2 = mean( chim ) ; % chi ˆ 2 averaged over a l l timeframes .

disp ( 'The inversion has finished' ) ;

i f lambda < lamlim | abs ( chi2 − 1) > e r r c h i
% Test i f shot r e c o n s t r u c t i o n i s acceptable .
disp ( ' ' ) ;
r e s u l t = [ 'Results: lambda=' num2str( lambda ) ', mean(chi2)=' num2str( chi2 ) ] ;
disp ( r e s u l t ) ;

disp ( [ 'WARNING: reconstruction seems to be wrong in shot ' . . .

int2str ( shot ) ' !' ] ) ;
disp ( ' ' ) ;
e r r f l a g = 1;

else

disp ( '  ' ) ;
r e s u l t = [ 'shot ' , int2str ( shot ) , ' OK. (lambda = ' , num2str( lambda ) , …
', chi2 = ' , num2str( chi2 ) , ')' ] ;
disp ( r e s u l t ) ;
disp ( ' ' ) ;
good = 1;

end

% SIXTH STEP : QUEST FOR THE X−POINT.
X = reshape ( Xinv , ny ∗ nx , length ( tvec ) ) ;
[K, J ] = find (X == max(max(X) ) ) ;
% K : Pixel containing the X−point .
Q = union ( 3 : 1 6 , 3 1 : 4 4 ) ;
% Q : P o t e n t i a l X−point r e g i o n .
K = i n t e r s e c t (K,Q) ;
% Limiting the X−point to the p o t e n t i a l X−point r e g i o n .
i f x b l i n d  % Optional X−point e x c l u s i o n

disp ( '  You have chosen to exlude the X-point' ) ;
disp ( '  Running a second inversion...' ) ;

disp ( '  ' ) ;
xpoint = 0;
i f isempty (K) % When the point with h i g h e s t e m i s s i v i t y l i e s a u t s i d e the

APPENDIX B. THE PROGRAM CODE 29

% p o t e n t i a l X−point r e g i o n .

disp ( 'Identification of X-point impossible' ) ;

xpoint = 1;

else

switch c u t o f f % Determination of the expanded X−point ( multiple p i x e l s )
case 0
[ p i x e l , J ] = find (X > ( l e v e l / 1 0 0 ) ∗ max(max(X) ) ) ;
% LEVEL s p e c i f i e s the i n f e r i o r r a d i a t i o n percentage l i m i t
% f o r p i x e l s p a r t i c i p a t i n g to the expanded X−point .
p i x e l = i n t e r s e c t ( p i x e l , Q) ;
case 1  % Overrides the automatic X−point determination .
p i x e l = [ 1 2 , 1 3 ] ;
case 2
p i x e l = [ 1 0 , 1 1 , 1 2 , 1 3 , 1 4 , 1 5 , 1 6 ] ;
case 3
p i x e l = [ 6 , 7 , 8 , 9 , 1 0 , 1 1 , 1 2 , 1 3 , 1 4 , 1 5 , 1 6 , 1 7 , 1 8 , 4 1 ] ;
case 4
p i x e l = [ 5 , 6 , 7 , 8 , 9 , 1 0 , 1 1 , 1 2 , 1 3 , 1 4 , 1 5 , 1 6 , 1 7 , 1 8 , 1 9 , …

2 0 , 4 0 , 4 1 , 4 2 , 4 3 ] ;

end

g = reshape ( Xinv / 1 0 0 , ny ∗ nx , length ( tvec ) ) ;
contX = T( inumdet , p i x e l ) ∗ g ( p i x e l , : ) ;
% Contribution of the expanded X−point to the chord
% b r i g h t n e s s Y.
Y = Y − contX ;
% New chord b r i g h t n e s s without e m i s s i v i t y from the expanded
% X−point .

APPENDIX B. THE PROGRAM CODE 30

[ Xinv , lambda , chim , inumdet ] = …
f a b r e c (Y, dY , xmesh , ymesh , tvec , inumdet , igoodchans , …
T , lambda init , l o o p f i s h e r , i t e r f i s h e r , i b o r d , l i m i t ) ;
% Second tomography i n v e r s i o n without expanded X−point c o n t r i b u t i o n s .
Xinv=Xinv ∗ 1 0 0 ;  % Conversion to [W/mˆ 3]
chi2=mean( chim ) ; % chi ˆ 2 averaged over a l l timeframes .
i f p l o t t i n g  % Shows a contour p l o t of the found e m i s s i v i t y .

figure ( 'Name' , 'Emissivity matrix for the last timeframe [W/m^3]' ) ;

contourf ( xmesh / 1 0 0 , ymesh / 1 0 0 , Xinv ( : , : , length ( tvec ) ) , 5 0 ) , …
shading f l a t , axis equal
colorbar ;

end

disp ( 'The second inversion has finished' ) ;

i f lambda < lamlim | abs ( chi2 − 1) > e r r c h i
% Test i f shot r e c o n s t r u c t i o n i s acceptable .
disp ( ' ' ) ;
r e s u l t = [ '  Results: lambda=' num2str( lambda ) ', mean(chi2)=' num2str( chi2 ) ] ;
disp ( r e s u l t ) ;

disp ( [ '  WARNING: second reconstruction not successful for shot ' int2str ( shot ) ' !' ] ) ;

disp ( ' ' ) ;
e r r f l a g =1;

else

disp ( '  ' ) ;
r e s u l t = [ 'shot ' , int2str ( shot ) , ' OK. (lambda = ' , num2str( lambda ) , …
', chi2 = ' , num2str( chi2 ) , ')' ] ;
good = 1;
disp ( r e s u l t ) ;
disp ( ' ' ) ;

end end

end

X = reshape ( Xinv , ny∗nx , length ( tvec ) ) ;
rm = xmesh / 1 0 0 ;  % Conversion to [m]

APPENDIX B. THE PROGRAM CODE 31

zm = ymesh / 1 0 0 ;  % Conversion to [m]
i f r a d i a t i o n  % Optional c a l c u l a t i o n of the r a d i a t i o n f r a c t i o n s i s s u e d i n s i d e and
% o u t s i d e of the LCFS s u r f a c e .

  disp ( [ 'Calculation of the LCFS fractions of the radiation started.' ] ) ;

[ Ptot , Pin , Pout , Pabove , Pbelow , tvec ] = …
fabpower (X, shot , tvec , xmesh , ymesh , 1 ) ;

end

mdsclose ;
i f saving  % Optinal saving of v a r i a b l e s on harddisk .
s a v e f i l e =['FABCAT_' , int2str ( shot ) ] ;

try

    eval ( [ 'save ' , [ savepath , s a v e f i l e ] , ' shot Xinv tvec FastPtot FastPbelow FastPxpoint good ' ] ) ;
    disp ( [ 'The following file has been created: ' , savepath s a v e f i l e ] ) ;

end end

% SEVENTH STEP : Fast c a l c u l a t i o n of the radiated power
for k = 1 : nx
Pinv ( : , k , : ) = Xinv ( : , k , : ) ∗ xmesh ( k ) ;
% Integrand f o r l a t e r i n t e g r a t i o n

end

f a c t =2∗pi∗dx∗dy/1 e6 ;  % I n t e g r a t i o n f a c t o r ( over t o r o i d a l r i n g ) .
P = reshape ( Pinv , nx∗ny , length ( tvec ) ) ;
% Pixel power matrix
FastPtot = sum(P)∗ f a c t ; % Sum over a l l p i x e l s
i f ˜ isempty (K)
FastPbelow = (sum(P( 1 : (K− 2 ) , : ) ) + sum(P( 2 9 : (K+ 2 6 ) , : ) ) ) ∗ f a c t ;
% Power radiated from the d i v e r t o r r e g i o n .

FastPxpoint = (sum(P(K−1:K+ 1 , : ) ) + sum(P(K+27:K+29,:))+P(K+ 5 6 , : ) ) ∗ f a c t ;
% Power radiated from the expanded X−point .

else

FastPbelow = NaN;
FastPxpoint = NaN;

end

i f MDS  % Optional saving to the MDS t r e e .
i f xpoint  % Comments f o r MDS node ’COMMENT’

       mtxt = [ 'FABCAT with X-point at pixel ' , num2str(K) ] ;

else

       mtxt = [ 'FABCAT without X-point at pixel ' , num2str(K) ] ;

end

i f isempty (K)  % No X−point found .

       mtxt = [ 'FABCAT, X-point determination failed' ] ;

end

f a b s t o r e ( shot , rm , zm , tvec , X, lambda , runsigma , chim , points , s e t d i f f ( 1 : 6 4 , inumdet ) , mtxt , r e s u l t , FastPtot , NaN, NaN, FastPxpoint , FastPbelow ) ;

end end

time end = cputime − debut ;
disp ( [ 'CPU runtime: ' num2str( time end ) , ' s.' ] ) ;
disp ( '  ' ) ;

B.1.2 fabsig.m

function [ framedata , frametbase , e r r f l a g , e r r s t r , goodchans ] = …
f a b s i g ( btshotnum , f cam , delchanuser , smoothmode , tmin , tmax , points , chantest )
% f u n c t i o n [ framedata , frametbase , e r r f l a g , e r r s t r , goodchans ] = …
%  f a b s i g ( btshotnum , f cam , delchanuser , smoothmode , tmin , tmax , points , chantest )

%

% Gets the bolometry s i g n a l s f o r the tomography a n a l y s i s , checks f o r bad

APPENDIX B. THE PROGRAM CODE 33

% bolometers and does data p r o c e s s i n g ( c a l i b r a t i o n ) . This i s the f a s t v e r s i o n .

%

% This f u n c t i o n has been r e w r i t t e n to s a t i s f y the needs of bolometry .
% I t i s based on btomodata .m written by Jan Mlynar in 1 9 9 9 .

%

% I f tmin i s a number , i t takes a l l time samples within tmin and tmax .
% Otherwise i t loads data only at LIUQE times

%

% History :

%

% C h r i s t i a n S c h l a t t e r , TPIV CRPP EPFL , February 2002
% major parts s t o l e n from the btomo package

%

% Release 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
i f nargin<5 tmin = NaN; end
i f nargin<7 points = NaN; end

fprintf ( 'FABCAT signal processing routine. Release 1.4 (2002/02/07)\n' ) ;
fprintf ( 'Loading of the bolometry signals.\n\n' ) ;

e r r f l a g = 0 ;  % Global e r r o r f l a g
e r r s t r = '' ;
s i g = [ ] ; tvec = [ ] ; delchans = [ ] ;
% Get d e f a u l t setup values
smoofac = 1 5 ;  % Bernard Joyce ’ s data treatment method i s
% f i t t i n g a polynomial over 2∗SMOOFAC timeframes .
% Default = 15 ( should not be changed ) .
f i l t f r e q = 5 0 ;  % Frequency of the butterworth 4 th order d i g i t a l
% f i l t e r used by C h r i s t i a n Deschenaux ’ s channel
% t e s t i n g and data treatment methods .
% Default = 5 0 [ Hz ] ( should not be changed ) .

%−%

% S i g n a l reading and s c a l i n g
i f isempty ( btshotnum )  % Test i f shotnumber was given .
e r r f l a g = 1;

   e r r s t r = ( [ 'You didn''t specify a shot number !!!' ] ) ;

mdsclose ;

return end APPENDIX B. THE PROGRAM CODE 34

mdsopen ( 'tcv_shot' , btshotnum ) ;  % Opening the MDS t r e e f o r reading .

b t a l l d a t a=mdsdata ( '\base::bolo:signals' ) ;

% Bolometry s i g n a l node .
tau = mdsdata ( '\base::bolo:tau' ) ; % Loading the bolometer thermal i n e r t i a
% constants ( f o r c a l i b r a t i o n ) .
i f isempty ( b t a l l d a t a )  % Check i f the shot e x i s t s in MDS.
e r r f l a g = 1;
e r r s t r = ( [ 'Sorry, no data available for shot ' , int2str ( btshotnum ) ] ) ;
mdsclose ;

return end

% ∗ ∗ ∗ I f TMIN was not s p e c i f i e d −> LIUQE times are loaded . ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
i f isnan ( tmin )

bttimebase=mdsdata ( [ 'dim_of(\base::bolo:signals)' ] ) ;

% Loading the NORMAL time base .

l i u q e t b a s e=mdsdata ( 'dim_of(\results::r_axis)' ) ;

% Loading the LIUQE time base .
mdsclose ;
i i =1;  % I n d i c a t o r f o r NORMAL times n e a r e s t to LIUQE times .
base = [ ] ;
tbase = [ ] ;
for k=1:length ( l i u q e t b a s e )
base ( i i )=max( find ( abs ( bttimebase−l i u q e t b a s e ( k))==min( abs ( bttimebase−l i u q e t b a s e ( k ) ) ) ) ) ;
% Index vector of NORMAL times n e a r e s t to LIUQE times .
i f i i ==1  % F i r s t n e a r e s t timeframe .
tbase ( i i )=bttimebase ( base ( i i ) ) ;
i i=i i +1;
else  % Following n e a r e s t timeframes .

APPENDIX B. THE PROGRAM CODE 35

i f bttimebase ( base ( i i ))> tbase ( i i −1)
tbase ( i i )=bttimebase ( base ( i i ) ) ;
i i=i i +1;

end end

end

frametbase=tbase ’ ;
[ vec , framebase , ib ]= i n t e r s e c t ( bttimebase , frametbase ) ;
% Times matching the time i n t e r v a l of the shots timeframes .
[ b t t o t s t e p s , bttotchan ]= size ( b t a l l d a t a ) ;
% ∗ ∗ ∗ I f TMIN was s p e c i f i e d −> NORMAL times are loaded . ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗

else

bttimebase=mdsdata ( [ 'dim_of(\base::bolo:signals)' ] ) ;

% Loading the NORMAL time base .
mdsclose ;
[ b t t o t s t e p s , bttotchan ]= size ( b t a l l d a t a ) ;
i f isnan ( tmax )  % I f TMAX was not s p e c i f i e d −> upper time l i m i t i s the
% l a t e s t timeframe f o r the given shot .
tmax = bttimebase ( length ( bttimebase )−1);

end

i f tmax > bttimebase ( length ( bttimebase ) )
tmax = bttimebase ( length ( bttimebase )−1);

end

i f ˜ isnan ( points )  % Reducing the number of timeframes to POINTS timeframes
% i f points have been submitted .
minstep = min( find ( bttimebase>=tmin ) ) ;

APPENDIX B. THE PROGRAM CODE 36

maxstep = max( find ( bttimebase<=tmax ) ) ;
framebase = round( minstep : ( maxstep−minstep )/( points −1): maxstep ) ;
else  % Taking a l l the timeframes matching the given time i n t e r v a l .
framebase = find ( bttimebase>=tmin & bttimebase<=tmax ) ;
points = length ( framebase ) ;

end

frametbase = bttimebase ( framebase ) ;

end %−%

% Channel t e s t i n g and d e l e t i n g
i f chantest ˜ = 0

fprintf ( 'Looking for bad channels:\n\n' ) ;

end

mdsopen ( 'tcv_shot' , btshotnum ) ;
b o l o s i g = mdsdata ( '\base::bolo:source' ) ;  % Reading of the unprocessed ( raw ) bolometry s i g n a l s
mdsclose ;
i f ˜ isempty ( delchanuser )  % Channels d e l e t e d by the user ( f l a g DELCHANUSER in
% ’ fabcat .m’

fprintf ( '  You have disabled the following channels:\n  ' ) ;

fprintf ( '%1.0f,' , delchanuser ’ ) ;
fprintf ( '\n\n' ) ;

end

% ∗ ∗ ∗ Automated d e l e t i n g of weird channels using Arno Refke ’ s method ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
% ∗ ∗ ∗ This i s the channel c o n s i s t e n c y check
i f ( ( chantest = = 2 ) | ( chantest == 3))
o f f s e t = mean( b o l o s i g ( 1 : 8 0 , : ) ) ;  % O f f s e t : Mean value of the s i g n a l s
b o l o s i g = b o l o s i g − repmat ( o f f s e t , b t t o t s t e p s , 1 ) ;
d i f b o l o = d i f f ( b o l o s i g ) ;  % Slope of the s i g n a l s .
gaga1 = find ( (max( b o l o s i g )−min( b o l o s i g ) ) < 0 . 5 ) ;  % S i g n a l value range check .
gaga2 = find (mean( b o l o s i g ) < 0 ) ;  % Negative s i g n a l check .

gaga3 = find ( (max( d i f b o l o )−min( d i f b o l o ) ) > 4 ) ;  % S i g n a l s l o p e check .
gaga4 = find ( abs ( o f f s e t ) > 9 ) ;  % S i g n a l o f f s e t check .
gaga = union ( gaga1 , gaga2 ) ;
gaga = union ( gaga , gaga3 ) ;
gaga = union ( gaga , gaga4 ) ;
i f ˜ isempty ( gaga )

  fprintf ( '  Arno Refke''s method disabled the following channels:\n  ' ) ;

fprintf ( '%1.0f,' , gaga ’ ) ;
fprintf ( '\n\n' ) ;

end

delchans=union ( gaga , delchanuser ) ;

end

% ∗ ∗ ∗ Automated d e l e t i n g of weird channels using C h r i s t i a n Deschenaux ’ s method ∗∗∗∗∗∗∗∗∗∗∗∗∗
% ∗ ∗ ∗ This i s the s i g n a l b a c k f i l t e r i n g method
i f ( ( chantest = = 1 ) | ( chantest == 3))
p = b t a l l d a t a ’ ;
amp=.01;
[ b , a]= butter ( 4 ,min( 1 , f i l t f r e q / ( 2 0 0 0 / 2 ) ) ) ;  % Design of a butterworth
% d i g i t a l f i l t e r of 4 th order .
[ py , px]= size (p ) ;
f l t r=zeros ( size ( bttotchan ) ) ;
for i =1: bttotchan
p( i , : ) = f i l t f i l t (b , a , p( i , : ) ) ;  % F i l t e r i n g of the s i g n a l s .
f l t r 1 ( i )=mean(p( i , [ 1 : min( 1 0 0 , b t t o t s t e p s ) ] ) ) ;
m=abs (mean(p( i , : ) ) ) ;
i f m>0  % F i l t e r e d s i g n a l mean value s i g n check .
f l t r 2 ( i )=mean( abs ( d i f f (p( i , : ) ) ) ) /m;  % F i l t e r e d s i g n a l s l o p e check

else

f l t r 2 ( i )=0;
end ;
end ;

APPENDIX B. THE PROGRAM CODE 38

f l t r 1=abs ( f l t r 1 )>0.1;
f l t r 2 =( f l t r 2 <(mean( f l t r 2 )∗amp) ) ;
f l t r=find ( f l t r 1 | f l t r 2 ) ;
i f ˜ isempty ( f l t r )

  fprintf ( '  Christian Deschenaux''s method disabled the following channels:\n  ' ) ;

fprintf ( '%1.0f,' , f l t r ’ ) ;
fprintf ( '\n\n' ) ;

end

delchans=union ( f l t r , delchans ) ;
delchans=union ( delchans , delchanuser ) ;

end

goodchans=s e t d i f f ( [ 1 : bttotchan ] , delchans ) ;

%−%

% Data treatment ( c a l i b r a t i o n )
i f smoothmode
sweep=smoofac ;

else

sweep =1;

end

i f framebase (1)< sweep+1  % Time c o m p a t i b i l i t y check ,
% to t e s t i f methods are a p p l i c a b l e .
e r r f l a g = 1;

e r r s t r = ( [ 'The smoothing requires higher start time.' ] ) ;

return

e l s e i f length ( framebase ) > b t t o t s t e p s − sweep
e r r f l a g = 1;

e r r s t r = ( [ 'The smoothing requires lower stop time.' ] ) ;

return end

rawdata=b t a l l d a t a ( framebase , goodchans ) ; % The ram , u n c a l i b r a t e d s i g n a l s .

framedata=NaN. ∗ ones ( size ( rawdata ) ) ;
% ∗ ∗ ∗ No data treatment ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
i f smoothmode == 0
framedata =1.78 e4∗rawdata ;  % C a l c u l a t i o n of the chord b r i g t h n e s s

end

% ∗ ∗ ∗ Bernard Joye ’ s method ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
% ∗ ∗ ∗ This i s the second order polynomial f i t t i n g method
i f smoothmode == 1

fprintf ( '  Smoothing using Bernard Joye''s method.\n' ) ;

o f f s e t = mean( b t a l l d a t a ( 1 : 8 0 , goodchans ) ) ;
basedata = 1.78 e4 ∗( b t a l l d a t a ( : , goodchans ) …
− repmat ( o f f s e t , b t t o t s t e p s , 1 ) ) ;
% Uncalibrated chord b r i g h t n e s s .
per acq =(bttimebase ( b t t o t s t e p s )−bttimebase ( 1 ) ) / ( b t t o t s t e p s −1);
% mean timestep time .
xvec=−smoofac : smoofac ;  % F i t t i n g index .
bolocor=NaN. ∗ ones ( 1 , b t t o t s t e p s ) ;
for stn =1:length ( framebase )  % F i t t i n g timeframe window .
k=framebase ( stn ) ;
yvec=basedata ( ( k−smoofac ) : ( k+smoofac ) , : ) ;
for chn=1:length ( goodchans )  % F i t t i n g s i g n a l s using a polinomial
% of second degree taking in account
% of the bolometers time constants TAU
p=polyfit ( xvec , yvec ( : , chn ) ’ , 2 ) ;
framedata ( stn , chn)=p ( 3 ) + …
tau ( goodchans ( chn ))∗ p (2)/ per acq ;

end end

end APPENDIX B. THE PROGRAM CODE 40

% ∗ ∗ ∗ C h r i s t i a n Deschenaux ’ s method ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
% ∗ ∗ ∗ This i s the b a c k f i l t e r i n g method
i f smoothmode == 2

   fprintf ( '  Smoothing using Christian Deschenaux''s method.\n' ) ;

[ b , a]= butter ( 4 ,min( 1 , f i l t f r e q / 1 0 0 0 ) ) ;  % Design of a butterworth
% d i g i t a l f i l t e r of 4 th order .
unity=ones ( size ( framebase ) ) ;
d i f f t=bttimebase ( framebase+unity)−bttimebase ( framebase−unity ) ;
% Neighboring fimeframes
% time i n t e r v a l s .
for ch=1:length ( goodchans )
p a l l= f i l t f i l t (b , a , b t a l l d a t a ( : , goodchans ( ch ) ) ) ; % F i l t e r i n g of the s i g n a l s .
p s e l=p a l l ( framebase ) ;
d i f f p=p a l l ( framebase+unity)− p a l l ( framebase−unity ) ;
% Derivation of the f i l t e r e d s i g n a l s .
framedata ( : , ch )=1.78 e4 ∗…
( p s e l+tau ( goodchans ( ch ))∗ d i f f p . / d i f f t ) ;
% C a l i b r a t i o n using time constants TAU.

end end

fprintf ( '\n' ) ;

B.1.3 fabrec.m

function [ g r e s , lamout , chi2out , a c t t ] = …
f a b r e c (Y, dY , xmesh , ymesh , tvec , a c t t , actdet , Tb , …
lambdai , ifishmax , iimax , i b o r d , l i m i t ) ;
% f u n c t i o n [ g r e s , lamout , chi2out , a c t t ] = …
%  f a b r e c (Y, dY , xmesh , ymesh , tvec , a c t t , actdet , Tb , …
%  lambdai , ifishmax , iimax , i b o r d , l i m i t ) ;

%

% Tomographic r e c o n s t r u c t i o n of the e m i s s i v i t y f o r the given
% s i g n a l s and geometry . The r e c o n s t r u c t i o n i s based on minimizing
% the Fisher information and uses a s i n g l e r e c o n s t r u c t i o n matrix
% which i s averaged over a l l timeframes .

%

% The r o u t i n e does a post−r e c o n s t r u c t i o n t e s t f o r bad channels and exclude
% them from the r e c o n s t r u c t i o n .

%

% Outputs : g r e s  : the e m i s s i v i t y array , s i z e ( ny x nx x length ( tvec ) )
%  lamout : the f i n a l smoothing f a c t o r LAMBDA
%  chi2out : The f i n a l chi ˆ 2 , length ( chi2out ) = lenght ( t s t e p s )
% a c t t : The array of enabled bolometer channels

%

% Please not that the average of chi2 = mean( chi2out ) i s targeted to 1 ,
% r e p r e s e n t i n g then the s o l u t i o n with maximum variance (dY ) .

%

% History :

%

% Jan Mlynarr , CRPP EPFL , May 2 0 0 1 ( f i r s t r e l e a s e )
% d i t o , January 2 0 0 2 ( c r o s s c h e c k of good channels )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )

fprintf ( 'FABCAT inversion routine. Release 1.4 (2002/02/07)\n' ) ;

while 1 % R e p e t i t i v e execution u n t i l no more weird channels are d i s c o v e r e d
% by the post−r e c o n s t r u c t i o n check ( see at the end of t h i s r o u t i n e )
lambda=lambdai ;
To=Tb( a c t t , : ) ;  % R e s t r i c t i o n s to enabled bolometers
sigmao=dY( a c t t ) ;
YY=Y( actdet , : ) ;
t s t e p s=length ( tvec ) ;
[ nl , npix ]= size (To ) ;
i f i b o r d  % no border c o n d i t i o n
c h i 2 t a r =( nl +1)∗ t s t e p s ;  % chi2 per a sample should be 1
% ( nl +1) one in plus − see the zero border below

else

c h i 2 t a r=nl ∗ t s t e p s ;  % chi2 per a sample should be 1

end

dx=xmesh(2)−xmesh ( 1 ) ;
dy=ymesh(2)−ymesh ( 1 ) ;
nx=length ( xmesh ) ;
ny=length ( ymesh ) ;
% START of the f a s t interim of of the cattcv package ’ makem 1 .m’
% only used when the no border c o n d i t i o n i s enabled .

e p s i l o n=sqrt ( eps ) ;
% Border p i x e l determination
i b o r d e r =[1: ny , ny+1:ny : ( ( ny ∗( nx −2))+1),(2:( nx−1))∗ny ,(1+( nx−1)∗ny ) : nx∗ny ] ;
nborder=length ( i b o r d e r ) ;
Ymax=max(YY) ;
D=max(YY. / repmat (Ymax, nl , 1 ) , e p s i l o n ∗ ones ( nl , t s t e p s ) ) ;
% Default model f o r the tought−f o r e m i s s i v i t y d i s t r i b u t i o n .
Tbar=sum(To ( : ) ) / ( nl+npix ) ;
Ep=D. / repmat (sum(To ’ ) ’ , 1 , t s t e p s ) ;
nn=find (sum(To)˜ =0);
Eca=(Ep’ ∗To ( : , nn ) ) . / repmat (sum(To ( : , nn ) ) , t s t e p s , 1 ) ;
% Totally f l a t model with zero border .
I0=sum( Eca ’ ) ;
Emean=I0 . / ( npix−nborder ) ;
g model=repmat (Emean , npix , 1 ) ;
g model ( i b o r d e r ,:)=1 e−9∗ones ( length ( i b o r d e r ) , t s t e p s ) ;
g model=g model . ∗ repmat (Ymax, npix , 1 ) ;
% END of f a s t interim of ” makem 1 .m”
tiny=1e −4;
dfadd=min( sigmao ) / 1 0 0 ; % R e l a t i v e e r r o r of the l i n e s of s i g h t
i p o c s =1;
rgmin=1e −2;
diam=eye ( npix ) ;  % diagonal
d i a r=diag ( ones ( 1 , npix−ny ) , ny ) ;  % to r e f e r e n c e r i g h t n e a r e s t neighbors
d i a l=diag ( ones ( 1 , npix−ny),−ny ) ;  % r e f e r e n c e l e f t n e a r e s t neighbors
diao=diag ( ones ( 1 , npix − 1 ) , 1 ) ;  % r e f e r e n c e upper n e a r e s t neighbors
diau=diag ( ones ( 1 , npix −1),−1);  % r e f e r e n c e lower n e a r e s t neighbors
% I n d i c e s f o r borders and c o i n s
i b l =2:ny−1;  % l e f t border
i b r =(2+(nx−1)∗ny ) : ( npix −1);  % r i g h t border
i b o =(2:( nx−1))∗ny ;  % upper border
i bu =1+(1:nx−2)∗ny ;  % lower border
i u l =1;  % lower l e f t corner
i o l=ny ;  % upper l e f t corner
i o r=npix ;  % upper r i g h t corner
i u r =1+(nx−1)∗ny ;  % lower r i g h t corner

i oben =[ i o l , i b o , i o r ] ;
i r e c h t s =[ i u r , i b r , i o r ] ;
% Evaluation of the g r a d i e n t s
Bx=−diam+d i a r ;  % f i n i t e−d i f f e r e n c e matrix r e p r e s e n t a t i o n s
By=−diam+diao ;  % of the p a r t i a l d e r i v a t e s .
By( i oben , : ) = diam ( i oben , : ) + diau ( i oben , : ) ;
Bx( i r e c h t s , : ) = diam ( i r e c h t s , : ) + d i a l ( i r e c h t s , : ) ;
Bx=Bx/dx ;
By=By/dy ;
% Normalization and d i v i s i o n by the variance
fo=Y( actdet , : ) ;
noemiss = [ ] ;  % border f i x e d to zero :
i f i b o r d  % border p i x e l s e m i s s i v i t y s e t to zero .
noemiss=find ( g model<=tiny ∗max(max( g model ) ) ) ;
Tadd=zeros ( 1 , npix ) ;
inx=find ( noemiss<npix ) ;
Tadd( noemiss ( inx ))= ones ( size ( inx ) ) ;
fadd=zeros ( 1 , t s t e p s ) ;
T=[To ; Tadd ] ;
f =[ fo ; fadd ] ;
sigma=[ sigmao ; dfadd ] ;

else

f = fo ;
T = To ;
sigma = sigmao ;

end

[ nl , npix ]= size (T) ;  % s i z e of T changes f o r the f i x e d border
fmax = max( f ) ;
f = f . / repmat ( fmax , size ( f , 1 ) , 1 ) ;  % Normalization
S1=diag ( 1 . / sigma ) ;  % D i v i s i o n by the variance
TT=zeros ( npix ) ;
Ts=S1∗T;
TT=Ts ’ ∗ Ts ;
f ( find ( f <0)) = zeros ( size ( find ( f <0))); % no n e g a t i v e s e m i s s i o n s :

f s = f . / repmat ( sigma , 1 , t s t e p s ) ;
lam=zeros ( ifishmax , iimax ) ;
chi2=lam ;
chi2det=zeros ( ifishmax , iimax , t s t e p s ) ;
chi2b=zeros ( ifishmax , t s t e p s ) ;
% ∗ ∗ ∗ Mainloop to minimize the Fisher information
for i =1: ifishmax
i f i ==1 % F i r s t loop
w = ones ( npix , 1 ) ;  % unity weight matrix
% The f i r s t i t e r a t i o n corresponds
% to l i n e a r r e g u l a r i z a t i o n
H = Bx’ ∗Bx+By’ ∗By ;  % The H−matrix
coeflam=trace (TT)/ trace (H) ;  % i n i t a l LAMBDA
lam (1,1)= lambda∗ coeflam ;
A=(TT+lam ( 1 , 1 ) ∗H ) ’ ;
Tpsinv=A\Ts ’ ;  % The M−matrix
g=Tpsinv∗ f s ;  % The e m i s s i v i t y
neg=find ( g <0);
g ( neg)=zeros ( size ( neg ) ) ;  % Negative e m i s s i v i t y i s s e t to zero
g ( noemiss)=zeros ( size ( noemiss ) ) ; % the same with noemission p i x e l s

else

gdummy=sum( g , 2 ) /max(sum( g , 2 ) ) ;
p e t i t=find (gdummy<rgmin ) ;
gdummy( p e t i t )=rgmin∗ ones ( size ( p e t i t ) ) ;
w=1./gdummy ;  % The weight matrix W(n ) , n>1
Ax=diag (w)∗Bx ;
Ay=diag (w)∗By ;
H=Bx’ ∗Ax+By’ ∗Ay ;  % H(n ) , n>1
lam ( i ,1)= lambdab ( i −1);
A=(TT+lam ( i , 1 ) ∗H ) ’ ;
Tpsinv=A\Ts ’ ;  % M(n ) , n>1
g=Tpsinv∗ f s ;
neg=find ( g <0);
g ( neg)=zeros ( size ( neg ) ) ;
g ( noemiss)=zeros ( size ( noemiss ) ) ;

end

lam ( i ,2)= lam ( i , 1 ) / 2 ;

f c=Ts∗g ;  % Reconstructed chord b r i g h t n e s s
dummy=f c−f s ;
chi2det ( i , 1 , : ) =sum(dummy. ∗dummy ) ’ ; % ( ifishmax ∗ iimax ∗ t s t e p s ) matrix
chi2 ( i ,1)=sum( chi2det ( i , 1 , : ) ) ;
% ∗ ∗ ∗ Subloop to f i n d LAMBDA which d e l i v e r s chi ˆ 2 −> t a r g e t chi ˆ 2
for i i =2: iimax
A=(TT+lam ( i , i i )∗H ) ’ ;
Tpsinv=A\Ts ’ ;
g=Tpsinv∗ f s ;
neg=find ( g <0);
g ( neg)=zeros ( size ( neg ) ) ;
g ( noemiss)=zeros ( size ( noemiss ) ) ;
f c=Ts∗g ;
dummy=f c−f s ;
chi2det ( i , i i , : ) =sum(dummy. ∗dummy ) ’ ;
chi2 ( i , i i )=sum( chi2det ( i , i i , : ) ) ;
i f i i=iimax break ; end  % I f l a s t subloop i t e r a t i o n reached
d e r i v =(lam ( i , i i )−lam ( i , i i −1))/( chi2 ( i , i i )−chi2 ( i , i i −1));
c o e f=lam ( i , i i )− d e r i v ∗ chi2 ( i , i i ) ;
lam ( i , i i +1)=d e r i v ∗ c h i 2 t a r+c o e f ; % The new lambda
i f lam ( i , i i +1)<=0 lam ( i , i i +1)=lam ( i , i i ) / 2 ; end
i i=i i +1;
end % of the subloop
[ noth , stepno ]=min( abs ( chi2 ( i ,:)− c h i 2 t a r ∗ ones ( size ( chi2 ( i , : ) ) ) ) ) ;
% f i n d lambda which makes chi2 c l o s e s t to c h i 2 t a r
chi2b ( i , : ) = chi2det ( i , stepno , : ) ;
lambdab ( i )=lam ( i , stepno ) ;  % Best LAMBDA so f a r
A=(TT+lambdab ( i )∗H ) ’ ;
Tpsinv=A\Ts ’ ;  % Best M−matrix so f a r
g=Tpsinv∗ f s ;  % Best e m i s s i v i t y r e c o n s t r u c t i o n so f a r
neg=find ( g <0);
g ( neg)=zeros ( size ( neg ) ) ;
g ( noemiss)=zeros ( size ( noemiss ) ) ;
end % of the main loop
lambda=lam/ coeflam ;  % ( ifishmax ∗ iimax ) matrix , d e s c r i b e s the i t e r a t i v e
% e v o l u t i o n of smoothing

APPENDIX B. THE PROGRAM CODE 46

lamout=lambda ( ifishmax , iimax ) ;  % Final LAMBDA
chi2out=chi2b ( ifishmax , : ) / nl ;  % Final chi ˆ 2 vector
g r e s=reshape ( g . ∗ repmat ( fmax , npix , 1 ) , ny , nx , length ( tvec ) ) ;
% Final e m i s s i v i t y
% Post−r e c o n s t r u c t i o n bad channel d e t e c t i o n
mdum=mean(dummy , 2 ) ;  % Mean chi
stdd=std (mdum) ;  % standard d e v i a t i o n of mean chi
goodch=find ( abs (mdum)< l i m i t ∗ stdd ) ;  % LIMIT determines the c o n f i d e n t i a l channels
i f length ( goodch)==length (mdum) break ; end
delch=find ( abs (mdum)>=l i m i t ∗ stdd ) ;
fprintf ( '\n' ) ;

fprintf ( '  Channels deleted by the chi2 test : ' ) ;

fprintf ( '%1.0f,' , a c t t ( delch ’ ) ) ;
fprintf ( '\n\n' ) ;
a c t t=a c t t ( goodch ) ;
actdet=actdet ( goodch ) ;

end

B.1.4 fabstore.m

function f a b s t o r e ( shot , rm , zm , time , X, lambda , Sigma , …
l o g c h i 2 , smoothing , nchords , mtxt , ctxt , Ptot , …
Pin , Pout , Pabove , Pbelow ) ;
% f u n c t i o n f a b s t o r e ( shot , rm , zm , time , X, lambda , Sigma , …
% l o g c h i 2 , smoothing , nchords , mtxt , ctxt , Ptot , …
%  Pin , Pout , Pabove , Pbelow ) ;

%

% Stores the r e s u l t s of the r e c o n s t r u c t i o n in the MDS t r e e .
% ATTENTION: whatever may be there w i l l be overwritten .

%

% This r o u t i n e has been taken from the b o l o t i package .

%

% Arguments :

%

%  shot  the shotnumber  [ 1 x 1 ]
%  rm  p i x e l c o o r d i n a t e s  [ 1 x nr ]
%  zm  idem  [ 1 x nz ]
%  time  times f o r the s l i c e s i n v e r t e d  [ 1 x timesteps ]
%  X  normalised e m i s s i v i t y ( max=1)  [ n p i x e l s x timesteps ]
%  lambda  smoothing f a c t o r of r e c o n s t r .  [ 1 x timesteps ]

APPENDIX B. THE PROGRAM CODE 47

%  Sigma  data e r r o r in %  [ 1 x 1 ]
%  l o g c h i 2  log10 ( chi2 / length ( nchords ) )  [ 3 x timesteps ]
%  smoothing  number of points averaged in  [ 1 x 1 ]
%  bolo−data treatment
%  nchords  number of bolo−chords used  [ 1 x 1 ]
%  mtxt  comment on i n v e r s i o n method  s t r i n g
%  Ptot  t o t a l radiated power [ Pin+Pout ] [ 1 x timesteps ]
%  Pin  radiated power i n s i d e LCFS  [ 1 x timesteps ]
%  Pout  radiated power o u t s i d e LCFS  [ 1 x timesteps ]
%  Pabove  radiated power above LCFS  [ 1 x timesteps ]
%  Pbelow  radiated power below LCFS  [ 1 x timesteps ]

%

% All r a d i a t i o n given in [W]

%

% History :

%

% Arno Refke , CRPP EPFL, 2 0 of November , 1 9 9 8 ( f i r s t r e l e a s e )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
i f nargin˜ =17

     error ( 'fabstore.m: sorry, incorrect number of input arguments' )

return

end

cmd1 = 'build_signal(build_with_units(f_float($1),"W"),*,f_float($2))' ;

mdsopen ( 'results' , shot )

mdsput ( '\results::btomo:method' , mtxt , 't' ) ;
mdsput ( '\results::btomo:comment' , ctxt , 't' ) ;
mdsput ( '\results::btomo:rmesh' , rm , 'f' ) ;
mdsput ( '\results::btomo:zmesh' , zm , 'f' ) ;
mdsput ( '\results::btomo:nchord' , nchords , 'f' ) ;
mdsput ( '\results::btomo:lambda' , lambda , 'f' ) ;
mdsput ( '\results::btomo:sigma' , Sigma , 'f' ) ;
mdsput ( '\results::btomo:chi_squared' , l o g c h i 2 , 'f' ) ;
mdsput ( '\results::btomo:smoothing' , smoothing , 'f' ) ;
mdsput ( '\results::btomo:time' , time , 'f' ) ;
mdsput ( '\results::btomo:emissivity' , . . . . .
         'Build_signal($1,*,*,\results::btomo:time)' , 'x' ,X) ;
mdsput ( '\results::btomo:prad_above' , cmd1 , 'x' , Pabove , time ) ;
mdsput ( '\results::btomo:prad_below' , cmd1 , 'x' , Pbelow , time ) ;
mdsput ( '\results::btomo:prad_in' , cmd1 , 'x' , Pin , time ) ;
mdsput ( '\results::btomo:prad_out' , cmd1 , 'x' , Pout , time ) ;

APPENDIX B. THE PROGRAM CODE 48

mdsput ( '\results::btomo:prad_tot' , cmd1 , 'x' , Ptot , time ) ;

mdsclose ;

return

B.2 Optional files

B.2.1 fabpower.m

function [ Ptot , Pin , Pout , Pabove , Pbelow , tvec ] = …
fabpower (X, shot , tvec , rm , zm , i p l o t ) ;
% f u n c t i o n [ Ptot , Pin , Pout , Pabove , Pbelow , tvec ] = …
%  fabpower (X, shot , tvec , rm , zm , i p l o t ) ;
%  ;
% C a l c u l a t e s the c e l l − and g r i d−points of a s e l e c t e d
% mesh ( rm , zm in [ cm ] ) f o r tomographic inverson of
% the bolometric data l y i n g i n s i d e / o u t s i d e the LCFS.
% Also determines the c u t p o i n t s of the l c f s with the tomogrid
% and c a l c u l a t e s Ptot , Pin , Pout , and Pdiv from the power
% radiated in−/o u t s i d e the l c f s from each c e l l

%

% This r o u t i n e was taken from the ’ btomo .m’ package

%

% History :

%

% Arno Refke , CRPP EPFL, 1 8 of September , 1 9 9 6 ( f i r s t r e l e a s e )
% Iwan J e r j e n , CRPP EPFL, 2 6 of August , 1 9 9 9 ( minor m o d i f i c a t i o n s )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
global h a f i g r a d  %handle of f i g u r e
i f nargin < 6
i p l o t = 0;
else i p l o t=i p l o t ;

end

% Normalize
c a l f=max(X ( : ) ) ; % s c a l a r
X=X/ c a l f ;
% Create vector Pm to every mesh−point
nx = length (rm ) ;  % number of x−c e l l s
ny = length (zm ) ;  % number of y c e l l s
dx = (rm(2)−rm ( 1 ) ) ;  % c e l l s i z e in x

dy = (zm(2)−zm ( 1 ) ) ;  % c e l l s i z e in y
xmax = rm( nx)+dx /2;  % maximum x mesh point
xmin = rm(1)−dx /2;  % minimum x mesh point
ymax = zm( ny)+dy /2;  % maximum y mesh point
ymin = ymax−ny∗dy ;  % minimum y mesh point
A c e l l = dx∗dy ;  % c e l l area ;
xgrid = [ xmin : dx : xmax +0.001];  % vector of x g r i d points
ygrid = f l i p l r ( [ ymin : dy : ymax + 0 . 0 0 1 ] ) ; % vector of y g r i d points
mx = length ( xgrid ) ;  % number of x meshgrid points
my = length ( ygrid ) ;  % number of y meshgrid−points
km = mx∗my;  % t o t a l number of mesh points
Pm=zeros (km , 3 ) ;
for j =1:my
Pm( (mx∗( j −1)+1):( j ∗mx) , 1 ) = xgrid ’ ;

end

for j =1:mx
Pm( j :mx: (mx∗(my−1)+ j ) , 2 ) = ygrid ’ ;

end

% Create vector Pc to every c e l l c e n t e r−point
% Only x c e l l and kc are used
x c e l l = rm ;  % vector of x c e l l c e n t e r points
kc = nx∗ny ;  % number of c e l l p o i n t s
i f length ( c a l f ) == 1, c a l f = c a l f ∗ ones ( length ( tvec ) , 1 ) ; end ;
% Loop over a l l times s p e c i f i e d in the shot
for l =1:length ( tvec ) ;
% Vector L to every point of the LCFS
[ con xt , con yt , c x t , c y t , times ]= f a b l c f s ( shot , tvec ( l ) ) ;
I = find (˜ isnan ( con xt ) & ˜ isnan ( con yt ) ) ;
con xt = con xt ( I ) ;
con yt = con yt ( I ) ;
l c f s x = con xt ;  % vector of l c f s x−points
l c f s y = con yt ;  % vector of l c f s y−points
y Xpt = min( l c f s y ) ; % y − coordinate of the X−point ;
y Tpt = max( l c f s y ) ; % y − coordinate of the top−point from LCFS;
i = length ( l c f s x ) ;  % number of l c f s−points
L=zeros ( i , 2 ) ;
L ( : , 1 ) = l c f s x ;

L ( : , 2 ) = l c f s y ;
% Matrix Rm with v e c t o r s from each mesh point to each l c f s−point
Rm = zeros ( i ,2∗km) ;
for j =1:2:(2∗km−1)
Rm( : , j ) = L( : , 1 ) −Pm( ( j +1)/2,1);
Rm( : , j +1) = L( : , 2 ) −Pm( ( j +1)/2,2);

end

% Determine f o r each mesh point the angle phi m from the p o l a r c o o r d i n a t e s of Rm
phi m = zeros ( i ,km) ;
phi m = atan2 (Rm( : , 2 : 2 : 2 ∗ km) ,Rm( : , 1 : 2 : ( 2 ∗ km)−1));
phi mn =phi m ; % phi mn element [− pi , pi ]
for j =1: i ,
for k=1:km,
i f phi m ( j , k)<0,
phi m ( j , k)=phi m ( j , k)+2∗pi ; % phi m element [ 0 , 2 ∗ pi ]

end end

end

dif phi m = max( phi m)−min( phi m ) ;
dif phi mn = max( phi mn)−min( phi mn ) ;
concavefactor =1.07;
for j =1:km
i f dif phi m ( j ) > concavefactor ∗pi & dif phi mn ( j ) > concavefactor ∗pi
Pm( j , 3 ) = 1 ;

else

Pm( j , 3 ) = 0 ;

end

end

[ Mout ] = find (Pm( : , 3 ) = = 0 ) ;
[ Min ] = find (Pm( : , 3 ) = = 1 ) ;
% Define d i f f e r e n t branches of the l c f s f o r use of i n t e r p o l a t i o n
[XMIN, N1]=min( con xt ) ;
[XMAX, N2]=max( con xt ) ;
[YMIN, N3]=min( con yt ) ;
[YMAX, N4]=max( con yt ) ;
l c f s x u = f l i p l r ( con xt (N2 : N1 ) ) ;
l c f s y u = f l i p l r ( con yt (N2 : N1 ) ) ;

l c f s x o = [ con xt (N1 : length ( con xt ) ) ; con xt ( 2 : N2 ) ] ;
l c f s y o = [ con yt (N1 : length ( con yt ) ) ; con yt ( 2 : N2 ) ] ;
i f N3>N4 ,
Nl = [N3 : length ( con yt ) 1 : N4 ] ;
Nr = [N4 : N3 ] ;
l c f s x l = con xt ( Nl ) ;
l c f s y l = con yt ( Nl ) ;
l c f s x r = con xt (Nr ) ;
l c f s y r = con yt (Nr ) ;
else  l c f s x l = con xt (N3 : N4 ) ;
l c f s y l = con yt (N3 : N4 ) ;
l c f s x r = f l i p l r ( [ con xt (N4 : length ( con xt ) ) ; con xt ( 2 : N3 ) ] ) ;
l c f s y r = f l i p l r ( [ con yt (N4 : length ( con yt ) ) ; con yt ( 2 : N3 ) ] ) ;

end

% Make sure that each branch i s monotonic
l c f s x o d i f = d i f f ( l c f s x o ) ;
[Mxo] = find ( l c f s x o d i f <0);
i f ˜ isempty (Mxo)
[ l c f s x o , I ] = sort ( l c f s x o ) ;
for m=1:length ( l c f s x o ) ,
l c f s y o (m) = l c f s y o ( I (m) ) ;

end end

l c f s x o d i f = d i f f ( l c f s x o ) ;
[ Nxo ] = find ( l c f s x o d i f ==0);
i f ˜ isempty (Nxo)
Nxo dif = d i f f (Nxo ) ;
[ nxo ] = find ( Nxo dif ==1);
i f length (Nxo)>0,
i f ˜ isempty ( nxo )
for m=1:length ( nxo ) ,
l c f s x o (Nxo( nxo (m) ) ) = ( l c f s x o (Nxo( nxo (m)))+ l c f s x o (Nxo( nxo (m))−1))/1.99999;

end end

for m=1:length (Nxo ) ,
i f Nxo(m) < length ( l c f s x o )−1,
l c f s x o (Nxo(m)+1) = ( l c f s x o (Nxo(m)+1)+ l c f s x o (Nxo(m)+2))/1.99999;
else l c f s x o (Nxo(m) ) = ( l c f s x o (Nxo(m))+ l c f s x o (Nxo(m)−1))/1.99999;

end end

end end

l c f s x u d i f = d i f f ( l c f s x u ) ;
[Mxu] = find ( l c f s x u d i f <0);
i f ˜ isempty (Mxu)
[ l c f s x u , I ] = sort ( l c f s x u ) ;

for m=1:length ( l c f s x u ) ,
l c f s y u (m) = l c f s y u ( I (m) ) ;

end end

l c f s x u d i f = d i f f ( l c f s x u ) ;
[ Nxu] = find ( l c f s x u d i f ==0);
i f ˜ isempty (Nxu)
Nxu dif = d i f f (Nxu ) ;
[ nxu ] = find ( Nxu dif ==1);
i f length (Nxu)>0,
i f ˜ isempty ( nxu )
for m=1:length ( nxu ) ,
l c f s x u (Nxu( nxu (m) ) ) = ( l c f s x u (Nxu( nxu (m)))+ l c f s x u (Nxu( nxu (m))−1))/1.99999;

end end

for m=1:length (Nxu ) ,
i f Nxu(m) < length ( l c f s x u )−1,
l c f s x u (Nxu(m)+1) = ( l c f s x u (Nxu(m)+1)+ l c f s x u (Nxu(m)+2))/1.99999;
else l c f s x u (Nxu(m) ) = ( l c f s x u (Nxu(m))+ l c f s x u (Nxu(m)−1))/1.99999;

end end

end end

l c f s y l d i f = d i f f ( l c f s y l ) ;
[ Myl ] = find ( l c f s y l d i f <0);
i f ˜ isempty (Myl)
[ l c f s y l , I ] = sort ( l c f s y l ) ;
for m=1:length ( l c f s y l ) ,
l c f s x l (m) = l c f s x l ( I (m) ) ;

end end

l c f s y l d i f = d i f f ( l c f s y l ) ;
[ Nyl ] = find ( l c f s y l d i f ==0);
i f ˜ isempty ( Nyl )
N y l d i f = d i f f ( Nyl ) ;
[ nyl ] = find ( N y l d i f ==1);
i f length ( Nyl )>0,
i f ˜ isempty ( nyl )
for m=1:length ( nyl ) ,
l c f s y l ( Nyl ( nyl (m) ) ) = ( l c f s y l ( Nyl ( nyl (m)))+ l c f s y l ( Nyl ( nyl (m))−1))/1.99999;

end end

for m=1:length ( Nyl ) ,
i f Nyl (m) < length ( l c f s y l )−1,
l c f s y l ( Nyl (m)+1) = ( l c f s y l ( Nyl (m)+1)+ l c f s y l ( Nyl (m)+2))/1.99999;
else l c f s y l ( Nyl (m) ) = ( l c f s y l ( Nyl (m))+ l c f s y l ( Nyl (m)−1))/1.99999;

end end

end end

l c f s y r d i f = d i f f ( l c f s y r ) ;
[ Myr] = find ( l c f s y r d i f <0);
i f ˜ isempty (Myr)
[ l c f s y r , I ] = sort ( l c f s y r ) ;
for m=1:length ( l c f s y r ) ,
l c f s x r (m) = l c f s x r ( I (m) ) ;

end end

l c f s y r d i f = d i f f ( l c f s y r ) ;
[ Nyr ] = find ( l c f s y r d i f ==0);
i f ˜ isempty ( Nyr )
Nyr dif = d i f f ( Nyr ) ;
[ nyr ] = find ( Nyr dif ==1);
i f length ( Nyr)>0,
i f ˜ isempty ( nyr )
for m=1:length ( nyr ) ,
l c f s y r ( Nyr ( nyr (m) ) ) = ( l c f s y r ( Nyr ( nyr (m)))+ l c f s y r ( Nyr ( nyr (m))−1))/1.99999;

end end

for m=1:length ( Nyr ) ,
i f Nyr (m) < length ( l c f s y r )−1,
l c f s y r ( Nyr (m)+1) = ( l c f s y r ( Nyr (m)+1)+ l c f s y r ( Nyr (m)+2))/1.99999;
else l c f s y r ( Nyr (m) ) = ( l c f s y r ( Nyr (m))+ l c f s y r ( Nyr (m)−1))/1.99999;

end end

end end

% C a l c u l a t i o n of Power radiated i n s i d e / o u t s i d e l c f s
pin ( l ) = 0 ;
pout ( l ) = 0 ;
pdiv ( l ) = 0 ;
ptop ( l ) = 0 ;
P in ( l ) = 0 ;
P out ( l ) = 0 ;
P div ( l ) = 0 ;
P top ( l ) = 0 ;
c e l l i n o u t = 2∗ ones ( kc ) ;
emiss = reshape (X ( : , l ) , length (zm ) , length (rm ) ) ;
emiss = flipud ( emiss ) ;
for j =1:ny
for k=1:nx
i f Pm( ( j −1)∗mx+k,3)==1 & Pm( ( j −1)∗mx+k+1,3)==1 & Pm( j ∗mx+k,3)==1 & Pm( j ∗mx+k+1,3)==1,
c e l l i n o u t ( ( j −1)∗nx+k ) = 1 ;

pin ( l ) = pin ( l ) + emiss ( j , k )∗ x c e l l ( k ) ;
e l s e i f Pm( ( j −1)∗mx+k,3)==0 & Pm( ( j −1)∗mx+k+1,3)==0 & Pm( j ∗mx+k,3)==0 & Pm( j ∗mx+k+1,3)==0,
c e l l i n o u t ( ( j −1)∗nx+k ) = 0 ;
pout ( l ) = pout ( l ) + emiss ( j , k )∗ x c e l l ( k ) ;
% Determines the c u t l i n e s of each c e l l with the l c f s
e l s e i f Pm( ( j −1)∗mx+k,3)==0 & Pm( ( j −1)∗mx+k+1,3)==0 & Pm( j ∗mx+k,3)==1 & Pm( j ∗mx+k+1,3)==1,
x cut (1) = xgrid ( k ) ;
y cut (1) = interp1 ( l c f s x o , l c f s y o , x cut ( 1 ) ) ;
x cut (2) = xgrid ( k+1);
y cut (2) = interp1 ( l c f s x o , l c f s y o , x cut ( 2 ) ) ;
Ain = ((min( y cut)−ygrid ( j +1))∗dx + abs ( d i f f ( y cut ))∗ dx / 2 ) / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Ain ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Ain ) ;

e l s e i f Pm( ( j −1)∗mx+k,3)==1 & Pm( ( j −1)∗mx+k+1,3)==1 & Pm( j ∗mx+k,3)==0 & Pm( j ∗mx+k+1,3)==0,
x cut (1) = xgrid ( k+1);
y cut (1) = interp1 ( l c f s x u , l c f s y u , x cut ( 1 ) ) ;
x cut (2) = xgrid ( k ) ;
y cut (2) = interp1 ( l c f s x u , l c f s y u , x cut ( 2 ) ) ;
Aout = ((min( y cut)−ygrid ( j +1))∗dx + abs ( d i f f ( y cut ))∗ dx / 2 ) / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Aout ) ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Aout ;
e l s e i f Pm( ( j −1)∗mx+k,3)==0 & Pm( ( j −1)∗mx+k+1,3)==1 & Pm( j ∗mx+k,3)==0 & Pm( j ∗mx+k+1,3)==1,
y cut (1) = ygrid ( j +1);
x cut (1) = interp1 ( l c f s y l , l c f s x l , y cut ( 1 ) ) ;
y cut (2) = ygrid ( j ) ;
x cut (2) = interp1 ( l c f s y l , l c f s x l , y cut ( 2 ) ) ;
Ain = ((min( x cut)−xgrid ( k ))∗ dy + abs ( d i f f ( x cut ))∗ dy / 2 ) / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Ain ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Ain ) ;
e l s e i f Pm( ( j −1)∗mx+k,3)==1 & Pm( ( j −1)∗mx+k+1,3)==0 & Pm( j ∗mx+k,3)==1 & Pm( j ∗mx+k+1,3)==0,
y cut (1) = ygrid ( j ) ;
x cut (1) = interp1 ( l c f s y r , l c f s x r , y cut ( 1 ) ) ;
y cut (2) = ygrid ( j +1);
x cut (2) = interp1 ( l c f s y r , l c f s x r , y cut ( 2 ) ) ;
Aout = ((min( x cut)−xgrid ( k ))∗ dy + abs ( d i f f ( x cut ))∗ dy / 2 ) / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Aout ) ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Aout ;
e l s e i f Pm( ( j −1)∗mx+k,3)==0 & Pm( ( j −1)∗mx+k+1,3)==0 & Pm( j ∗mx+k,3)==0 & Pm( j ∗mx+k+1,3)==1,
y cut (1) = ygrid ( j +1);
x cut (1) = interp1 ( l c f s y l , l c f s x l , y cut ( 1 ) ) ;
x cut (2) = xgrid ( k+1);

y cut (2) = interp1 ( l c f s x o , l c f s y o , x cut ( 2 ) ) ;
Ain = abs ( d i f f ( x cut ))∗ abs ( d i f f ( y cut ) ) / 2 / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Ain ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Ain ) ;
e l s e i f Pm( ( j −1)∗mx+k,3)==0 & Pm( ( j −1)∗mx+k+1,3)==0 & Pm( j ∗mx+k,3)==1 & Pm( j ∗mx+k+1,3)==0,
y cut (2) = ygrid ( j +1);
x cut (2) = interp1 ( l c f s y r , l c f s x r , y cut ( 2 ) ) ;
x cut (1) = xgrid ( k ) ;
y cut (1) = interp1 ( l c f s x o , l c f s y o , x cut ( 1 ) ) ;
Ain = abs ( d i f f ( x cut ))∗ abs ( d i f f ( y cut ) ) / 2 / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Ain ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Ain ) ;
e l s e i f Pm( ( j −1)∗mx+k,3)==1 & Pm( ( j −1)∗mx+k+1,3)==0 & Pm( j ∗mx+k,3)==0 & Pm( j ∗mx+k+1,3)==0,
y cut (1) = ygrid ( j ) ;
x cut (1) = interp1 ( l c f s y r , l c f s x r , y cut ( 1 ) ) ;
x cut (2) = xgrid ( k ) ;
y cut (2) = interp1 ( l c f s x u , l c f s y u , x cut ( 2 ) ) ;
Ain = abs ( d i f f ( x cut ))∗ abs ( d i f f ( y cut ) ) / 2 / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Ain ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Ain ) ;
e l s e i f Pm( ( j −1)∗mx+k,3)==0 & Pm( ( j −1)∗mx+k+1,3)==1 & Pm( j ∗mx+k,3)==0 & Pm( j ∗mx+k+1,3)==0,
y cut (2) = ygrid ( j ) ;
x cut (2) = interp1 ( l c f s y l , l c f s x l , y cut ( 2 ) ) ;
x cut (1) = xgrid ( k+1);
y cut (1) = interp1 ( l c f s x u , l c f s y u , x cut ( 1 ) ) ;
Ain = abs ( d i f f ( x cut ))∗ abs ( d i f f ( y cut ) ) / 2 / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Ain ;
P out ( l ) = P out ( l )+ emiss ( j , k )∗ x c e l l ( k)∗(1−Ain ) ;
e l s e i f Pm( ( j −1)∗mx+k,3)==0 & Pm( ( j −1)∗mx+k+1,3)==1 & Pm( j ∗mx+k,3)==1 & Pm( j ∗mx+k+1,3)==1,
x cut (1) = xgrid ( k ) ;
y cut (1) = interp1 ( l c f s x o , l c f s y o , x cut ( 1 ) ) ;
y cut (2) = ygrid ( j ) ;
x cut (2) = interp1 ( l c f s y l , l c f s x l , y cut ( 2 ) ) ;
Aout = abs ( d i f f ( x cut ))∗ abs ( d i f f ( y cut ) ) / 2 / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Aout ) ;
P out ( l ) = P out ( l )+ emiss ( j , k )∗ x c e l l ( k )∗ Aout ;
e l s e i f Pm( ( j −1)∗mx+k,3)==1 & Pm( ( j −1)∗mx+k+1,3)==0 & Pm( j ∗mx+k,3)==1 & Pm( j ∗mx+k+1,3)==1,
x cut (2) = xgrid ( k+1);
y cut (2) = interp1 ( l c f s x o , l c f s y o , x cut ( 2 ) ) ;
y cut (1) = ygrid ( j ) ;
x cut (1) = interp1 ( l c f s y r , l c f s x r , y cut ( 1 ) ) ;
Aout = abs ( d i f f ( x cut ))∗ abs ( d i f f ( y cut ) ) / 2 / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Aout ) ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Aout ;

e l s e i f Pm( ( j −1)∗mx+k,3)==1 & Pm( ( j −1)∗mx+k+1,3)==1 & Pm( j ∗mx+k,3)==1 & Pm( j ∗mx+k+1,3)==0,
x cut (1) = xgrid ( k+1);
y cut (1) = interp1 ( l c f s x u , l c f s y u , x cut ( 1 ) ) ;
y cut (2) = ygrid ( j +1);
x cut (2) = interp1 ( l c f s y r , l c f s x r , y cut ( 2 ) ) ;
Aout = abs ( d i f f ( x cut ))∗ abs ( d i f f ( y cut ) ) / 2 / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Aout ) ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Aout ;
e l s e i f Pm( ( j −1)∗mx+k,3)==1 & Pm( ( j −1)∗mx+k+1,3)==1 & Pm( j ∗mx+k,3)==0 & Pm( j ∗mx+k+1,3)==1,
x cut (2) = xgrid ( k ) ;
y cut (2) = interp1 ( l c f s x u , l c f s y u , x cut ( 2 ) ) ;
y cut (1) = ygrid ( j +1);
x cut (1) = interp1 ( l c f s y l , l c f s x l , y cut ( 1 ) ) ;
Aout = abs ( d i f f ( x cut ))∗ abs ( d i f f ( y cut ) ) / 2 / A c e l l ;
P in ( l ) = P in ( l ) + emiss ( j , k )∗ x c e l l ( k)∗(1−Aout ) ;
P out ( l ) = P out ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Aout ;

end end

end

% C a l c u l a t i o n of power radiated in the d i v e r t o r ( below the X−point )
for j =1:ny
for k=1:nx
i f Pm( ( j −1)∗mx+k,2)< y Xpt & Pm( ( j −1)∗mx+k+1,2)<y Xpt & Pm( j ∗mx+k,2)< y Xpt & Pm( j ∗mx+k+1,2)<y Xpt ,
pdiv ( l ) = pdiv ( l ) + emiss ( j , k )∗ x c e l l ( k ) ;
e l s e i f Pm( ( j −1)∗mx+k,2)> y Xpt & Pm( ( j −1)∗mx+k+1,2)>y Xpt & Pm( j ∗mx+k,2)< y Xpt & Pm( j ∗mx+k+1,2)<y Xpt ,
x cut (1) = xgrid ( k ) ;
x cut (2) = xgrid ( k+1);
Adiv = d i f f ( x cut )∗( y Xpt − ygrid ( j +1))/ A c e l l ;
P div ( l ) = P div ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Adiv ;

end end

end

% C a l c u l a t i o n of power radiated above the LCFS
for j =1:ny
for k=1:nx
i f Pm( ( j −1)∗mx+k,2)> y Tpt & Pm( ( j −1)∗mx+k+1,2)>y Tpt & Pm( j ∗mx+k,2)> y Tpt & Pm( j ∗mx+k+1,2)>y Tpt ,
ptop ( l ) = ptop ( l ) + emiss ( j , k )∗ x c e l l ( k ) ;
e l s e i f Pm( ( j −1)∗mx+k,2)> y Tpt & Pm( ( j −1)∗mx+k+1,2)>y Tpt & Pm( j ∗mx+k,2)< y Tpt & Pm( j ∗mx+k+1,2)<y Tpt ,
x cut (1) = xgrid ( k ) ;

x cut (2) = xgrid ( k+1);
Atop = d i f f ( x cut )∗( ygrid ( j )−y Tpt )/ A c e l l ;
P top ( l ) = P top ( l ) + emiss ( j , k )∗ x c e l l ( k )∗ Atop ;

end end

end

[ IN ] = find ( c e l l i n o u t == 1);
[OUT] = find ( c e l l i n o u t == 0);
% Power radiated i n s i d e / o u t s i d e LCFS from c e l l s
P in ( l ) = c a l f ( l )∗ P in ( l )∗2∗ pi∗dx∗dy ; % p a r t i a l l y i n s i d e l c f s
P out ( l ) = c a l f ( l )∗ P out ( l )∗2∗ pi∗dx∗dy ; % p a r t i a l l y o u t s i d e l c f s
P div ( l ) = c a l f ( l )∗ P div ( l )∗2∗ pi∗dx∗dy ; % p a r t i a l l y below  l c f s
P top ( l ) = c a l f ( l )∗ P top ( l )∗2∗ pi∗dx∗dy ; % p a r t i a l l y above  l c f s
pin ( l ) = c a l f ( l )∗ pin ( l )∗2∗ pi∗dx∗dy ; % t o t a l l y i n s i d e l c f s
pout ( l ) = c a l f ( l )∗ pout ( l )∗2∗ pi∗dx∗dy ; % t o t a l l y o u t s i d e l c f s
pdiv ( l ) = c a l f ( l )∗ pdiv ( l )∗2∗ pi∗dx∗dy ; % t o t a l l y below  l c f s
ptop ( l ) = c a l f ( l )∗ ptop ( l )∗2∗ pi∗dx∗dy ; % t o t a l l y above  l c f s
Pin ( l ) = ( pin ( l ) + P in ( l ))∗1 e −6;  % i n s i d e l c f s rad . power
Pout ( l ) = ( pout ( l )+ P out ( l ))∗1 e −6;  % o u t s i d e l c f s rad . power
Pdiv ( l ) = ( pdiv ( l )+ P div ( l ))∗1 e −6;  % below  l c f s rad . power
Ptop ( l ) = ( ptop ( l )+ P top ( l ))∗1 e −6;  % above  l c f s rad . power
Ptot ( l ) = ( Pin ( l ) + Pout ( l ) ) ;  % t o t a l l y radiated power

end

Pabove = Ptop ;
Pbelow = Pdiv ;
% Plots radiated power f r a c t i o n s ( Ptot , Pin , Pout and Pdiv )
% Calculated from bolometric i n v e r s i o n r o u t i n e
i f i p l o t

  h a f i g r a d=figure ( 'Name' , 'btomorun radiation fractions window' , . . .
         'NumberTitle' , 'off' ) ;

c l f

hmenu=uimenu( 'Parent' , h a f i g r a d , …

           'Label' , 'E&XIT' , . . .
           'Callback' , 'close(gcf)' ) ;

plot ( tvec , Ptot /1 e +3,'b' )
hold on
plot ( tvec , Pin /1 e +3,'r' )
plot ( tvec , Pout/1 e +3,'g' )
plot ( tvec , Pdiv /1 e +3,'-.m' )
plot ( tvec , Ptop/1 e +3,'--c' )
t i t l e ( [ '# ' , int2str ( shot ) , ' bolometric inversion' ] ) ;

  xlabel ( 'time  [s]' ) ;

APPENDIX B. THE PROGRAM CODE 58

   ylabel ( 'Prad  [kW]' ) ;
   legend ( 'Ptot' , 'Pin' , 'Pout' , 'Pbelow' , 'Pabove' , 0 ) ;

end

B.2.2 fablcfs.m

function [ con xt , con yt , c x t , c y t , c times ] = …
f a b l c f s ( shot , times )
% f u n c t i o n [ con xt , con yt , c x t , c y t , c times ] = …
%  f a b l c f s ( shot , times )

%

% This r o u t i n e r e t u r n s matrices of the contours of the
% LCFS from LIUQE f o r #shot as w e l l as the magnetic a x i s .
% ’ c times ’ contains the LIUQE times which were n e a r e s t
% to the values s p e c i f i e d in the vector ” times ” .

%

% s i z e of con xt i s [ npts contour x length ( c times ) ]
% s i z e of c y t i s [ 1 x length ( c times ) ] .

%

% This r o u t i n e has been taken form the ’ btomo .m’ package

%

% History :

%

% Matthias Anton , CRPP EPFL, 1 9 9 5 ( f i r s t r e l e a s e )
% Iwan J e r j e n , CRPP EPFL, 2 3 of August , 1 9 9 9 ( minor m o d i f i c a t i o n s )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
sshot=int2str ( shot ) ;
% GET THE DATA FROM MDS
% check a v a i l a b i l i t y of data
shot=mdsopen ( 'tcv_shot' , shot ) ;

pts=mdsdata ( '\results::npts_contour')−1;

% check whether data i s a v a i l a b l e f o r t h i s shot
i f ( pts (1) == −1)

   disp ( [ 'NO LIUQE DATA AVAILABLE FOR SHOT: ' , sshot ] ) ;

return ;

end

% get v a l i d timebase

APPENDIX B. THE PROGRAM CODE 59

tbase=mdsdata ( 'dim_of(\results::r_axis)' ) ;

% ensure that ’ time ’ i s within range
i f min( times)<min( tbase )
c time1 =1.001∗min( tbase ) ; s c t i m e s 1=num2str( c time1 ) ;

  disp ( [ 'WARNING: min time   out of range, using T =' s c t i m e s 1 '[s] instead' ] ) ;

end

i f max( times)>max( tbase )
c time2 =0.999∗max( tbase ) ; s c t i m e s 2=num2str( c time2 ) ;

  disp ( [ 'WARNING: max time out of range, using T =' s c t i m e s 2 '[s] instead' ] ) ;

end

i i =1;
for k=1:length ( times )
indext=max( find ( abs ( times ( k)−tbase)==min( abs ( times ( k)−tbase ) ) ) ) ;
i f i i ==1
c times ( i i )=tbase ( indext ) ;
i i=i i +1;

else

i f tbase ( indext)> c times ( i i −1)
c times ( i i )=tbase ( indext ) ;
i i=i i +1;

end end

end

con xt=mdsdata ( '\results::r_contour[*,$]' , c times ) ;
con yt=mdsdata ( '\results::z_contour[*,$]' , c times ) ;

con xt=con xt ∗100;
con yt=con yt ∗100;
c x t =(mdsdata ( '\results::r_axis[$]' , c times ) ∗ 1 0 0 ) ;

APPENDIX B. THE PROGRAM CODE 60

c y t=mdsdata ( '\results::z_axis[$]' , c times )∗100;
mdsclose ;
i f length ( c times )>1
con xt ( find ( con xt==0))=NaN∗ ones ( size ( find ( con xt ==0)));
con yt ( find ( con xt==0))=NaN∗ ones ( size ( find ( con xt ==0)));

else

con xt=con xt ( find ( con xt ˜ =0));
con yt=con yt ( find ( con xt ˜ =0));

end return

B.2.3 fabsetetendue.m

function fabsetetendue
% f u n c t i o n fabsetetendue

%

% C a l c u l a t e s d e t e c t o r etendue f o r bolometer tomography .

%

% History :

%

% C h r i s t i a n S c h l a t t e r , TPIV CRPP EPFL , February 2 0 0 2 ( f i r s t r e l e a s e )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
disp ( ' ' ) ;

disp ( 'Etendue calculator for bolometry.' ) ;
disp ( 'Final release 1.4' ) ;
disp ( 'Latest modification on 2002/02/07.' ) ;

disp ( ' ' ) ;
cw=1;  % d e t e c t o r numbers cw = 1 : c l o c k w i s e
% cw = 0 : c o u n t e r c l o c k w i s e
fans = [ 1 1 1 1 1 1 1 1 ] ;  % camera switch
vangle = [ 9 0 0 0 0 0 0 0 − 9 0 ] ;  % angle of d e t e c t o r s u r f a c e normal
xpos = [ 8 8 1 2 3 . 5 1 2 3 . 5 1 2 3 . 5 1 2 3 . 5 1 2 3 . 5 1 2 3 . 5 8 8 ] ;
% x p o s i t i o n of the diaphragmas in [ cm ]
ypos =[81.5 45.5 45.5 −.25 −.25 −46 −46 −81.5];

% y p o s i t i o n of the diaphragmas in [ cm ]
% P o s i t i o n s of the d e t e c t o r s
rdet = [ 0 . 8 6 2 , 0 . 8 6 7 , 0 . 8 7 2 1 , 0 . 8 7 7 1 , 0 . 8 8 2 9 , 0 . 8 8 7 9 , 0 . 8 9 3 , 0 . 8 9 8 , …

1 . 2 7 8 , 1 . 2 8 0 6 , 1 . 2 8 3 3 , 1 . 2 8 6 1 , 1 . 2 8 8 3 , 1 . 2 8 9 3 , 1 . 2 9 0 2 , 1 . 2 9 1 2 ] ;

rdet (17:24)= f l i p l r ( rdet ( 9 : 1 6 ) ) ;
rdet (25:40)= rdet ( 9 : 2 4 ) ;
rdet (41:56)= rdet ( 9 : 2 4 ) ;
rdet (57:64)= f l i p l r ( rdet ( 1 : 8 ) ) ;
zdet = [ 0 . 9 0 5 7 , 0 . 9 0 6 3 , 0 . 9 0 6 9 , 0 . 9 0 7 5 ] ;
zdet (5:8)= f l i p l r ( zdet ( 1 : 4 ) ) ;
zdet ( 9 : 2 4 ) = [ 0 . 4 9 1 6 , 0 . 4 8 7 3 , 0 . 4 8 3 , 0 . 4 7 8 7 , 0 . 4 7 3 , 0 . 4 6 8 , 0 . 4 6 3 1 , 0 . 4 5 8 1 , …

0 . 4 5 1 9 , 0 . 4 4 6 9 , 0 . 4 4 2 , 0 . 4 3 7 , 0 . 4 3 1 3 , 0 . 4 2 7 , 0 . 4 2 2 7 , 0 . 4 1 8 4 ] ;

zdet (25:40)= zdet ( 9 : 2 4 ) − 0 . 4 5 7 5 ;
zdet (41:56)= zdet ( 9 : 2 4 ) − 0 . 9 1 5 ;
zdet (57:64)=− zdet ( 1 : 8 ) ;
xdet=rdet ∗ 1 0 0 ; ydet=zdet ∗100;
d1 =1.5/10;  % d e t e c t o r width in cm
d2=4/10;  % d e t e c t o r length in cm
b1 = [ 2 . 6 2 . 2 2 . 2 2 . 2 2 . 2 2 . 2 2 . 2 2 . 6 ] / 5 ;
% aperture width in cm ( p o l o i d a l )
b2 = [ 1 0 8 8 8 8 8 8 1 0 ] / 5 ;  % aperture length in cm ( t o r o i d a l )
% C a l c u l a t i o n of aperture and d e t e c t o r p o s i t i o n s
nact=sum( fans ) ;
i a c t=find ( fans ) ;
ndet =8;
ncam=8;
xap=ones ( ndet , 1 ) ∗ xpos ;
xap=xap ( : ) ’ ;
yap=ones ( ndet , 1 ) ∗ ypos ;
yap=yap ( : ) ’ ;
be1=ones ( ndet , 1 ) ∗ b1 ;
be1=be1 ( : ) ’ ;
be2=ones ( ndet , 1 ) ∗ b2 ;
be2=be2 ( : ) ’ ;
x0=xdet−xap ;
y0=ydet−yap ;
i v e r t = [ 1 : 8 , 5 7 : 6 4 ] ’ ;

x0 ( i v e r t )=ydet ( i v e r t )−yap ( i v e r t ) ;
y0 ( i v e r t )=xdet ( i v e r t )−xap ( i v e r t ) ;
alpha=atan ( y0 . / x0 ) ;
rsquare=x0.ˆ 2+y0 . ˆ 2 ;
etend=d1∗d2∗be1 . ∗ be2 . ∗ cos ( alpha ) . ˆ 2 . / rsquare ;
% The bolometers etendue
save fabetendue . mat etend ; % Save to disk

disp ( 'File fabetendue.mat successfully created.' ) ;

disp ( ' ' ) ;

B.2.4 fabsetchord.m

function c a t s e t c h o r d
% f u n c t i o n c a t s e t c h o r d

%

% Detector viewing l i n e s c o o r d i n a t e s .

%

% Release 1 . 4 ( 0 7 / 0 2 / 2 0 0 2 )
disp ( ' ' ) ;

disp ( 'FABCAT Detector viewing line calculator for bolometry.' ) ;
disp ( 'Final Release 1.4.' ) ;
disp ( 'Latest modification on 2002/02/07.' ) ;

disp ( ' ' ) ;
xchord ( 1 , 0 1 : 0 9 ) = [  0 . 8 6 2 0  0 . 8 6 7 0  0 . 8 7 2 1  0 . 8 7 7 1  0 . 8 8 2 9  0 . 8 8 7 9  0 . 8 9 3 0  0 . 8 9 8 0  1 . 2 7 8 0 ] ;
xchord ( 2 , 0 1 : 0 9 ) = [  1 . 1 3 7 0  1 . 0 8 4 3  1 . 0 1 1 1  0 . 9 2 9 1  0 . 8 3 0 9  0 . 7 4 5 5  0 . 6 5 7 2  0 . 6 2 0 0  0 . 6 2 0 0 ] ;
xchord ( 1 , 1 0 : 1 8 ) = [  1 . 2 8 0 6  1 . 2 8 3 3  1 . 2 8 6 1  1 . 2 8 8 3  1 . 2 8 9 3  1 . 2 9 0 2  1 . 2 9 1 2  1 . 2 9 1 2  1 . 2 9 0 2 ] ;
xchord ( 2 , 1 0 : 1 8 ) = [  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0 ] ;
xchord ( 1 , 1 9 : 2 7 ) = [  1 . 2 8 9 3  1 . 2 8 8 3  1 . 2 8 6 1  1 . 2 8 3 3  1 . 2 8 0 6  1 . 2 7 8 0  1 . 2 7 8 0  1 . 2 8 0 6  1 . 2 8 3 3 ] ;
xchord ( 2 , 1 9 : 2 7 ) = [  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 7 2 6 1  0 . 8 1 8 5  0 . 8 8 8 4  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0 ] ;
xchord ( 1 , 2 8 : 3 6 ) = [  1 . 2 8 6 1  1 . 2 8 8 3  1 . 2 8 9 3  1 . 2 9 0 2  1 . 2 9 1 2  1 . 2 9 1 2  1 . 2 9 0 2  1 . 2 8 9 3  1 . 2 8 8 3 ] ;
xchord ( 2 , 2 8 : 3 6 ) = [  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0 ] ;
xchord ( 1 , 3 7 : 4 5 ) = [  1 . 2 8 6 1  1 . 2 8 3 3  1 . 2 8 0 6  1 . 2 7 8 0  1 . 2 7 8 0  1 . 2 8 0 6  1 . 2 8 3 3  1 . 2 8 6 1  1 . 2 8 8 3 ] ;
xchord ( 2 , 3 7 : 4 5 ) = [  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 8 9 4 3  0 . 8 2 5 6  0 . 7 3 4 7  0 . 6 2 0 0  0 . 6 2 0 0 ] ;
xchord ( 1 , 4 6 : 5 4 ) = [  1 . 2 8 9 3  1 . 2 9 0 2  1 . 2 9 1 2  1 . 2 9 1 2  1 . 2 9 0 2  1 . 2 8 9 3  1 . 2 8 8 3  1 . 2 8 6 1  1 . 2 8 3 3 ] ;
xchord ( 2 , 4 6 : 5 4 ) = [  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0 ] ;
xchord ( 1 , 5 5 : 6 3 ) = [  1 . 2 8 0 6  1 . 2 7 8 0  0 . 8 9 8 0  0 . 8 9 3 0  0 . 8 8 7 9  0 . 8 8 2 9  0 . 8 7 7 1  0 . 8 7 2 1  0 . 8 6 7 0 ] ;
xchord ( 2 , 5 5 : 6 3 ) = [  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 2 0 0  0 . 6 5 7 2  0 . 7 4 5 5  0 . 8 3 0 9  0 . 9 2 9 1  1 . 0 1 1 1  1 . 0 8 4 3 ] ;
xchord ( 1 , 6 4 : 6 4 ) = [  0 . 8 6 2 0 ] ;
xchord ( 2 , 6 4 : 6 4 ) = [  1 . 1 3 7 0 ] ;

ychord ( 1 , 0 1 : 0 9 ) = [  0 . 9 0 5 7  0 . 9 0 6 3  0 . 9 0 6 9  0 . 9 0 7 5  0 . 9 0 7 5  0 . 9 0 6 9  0 . 9 0 6 3  0 . 9 0 5 7  0 . 4 9 1 6 ] ;
ychord ( 2 , 0 1 : 0 9 ) = [ −0.4800 −0.6200 −0.7100 −0.7500 −0.7500 −0.7500 −0.7500 −0.4951 −0.0685];
ychord ( 1 , 1 0 : 1 8 ) = [  0 . 4 8 7 3  0 . 4 8 3 0  0 . 4 7 8 7  0 . 4 7 3 0  0 . 4 6 8 0  0 . 4 6 3 1  0 . 4 5 8 1  0 . 4 5 1 9  0 . 4 4 6 9 ] ;
ychord ( 2 , 1 0 : 1 8 ) = [  0 . 0 1 9 4  0 . 0 9 8 5  0 . 1 6 9 8  0 . 2 4 7 3  0 . 3 0 7 8  0 . 3 6 4 8  0 . 4 2 1 1  0 . 4 8 8 9  0 . 5 4 5 2 ] ;
ychord ( 1 , 1 9 : 2 7 ) = [  0 . 4 4 2 0  0 . 4 3 7 0  0 . 4 3 1 3  0 . 4 2 7 0  0 . 4 2 2 7  0 . 4 1 8 4  0 . 0 3 4 1  0 . 0 2 9 8  0 . 0 2 5 5 ] ;
ychord ( 2 , 1 9 : 2 7 ) = [  0 . 6 0 2 2  0 . 6 6 2 7  0 . 7 4 0 2  0 . 7 5 0 0  0 . 7 5 0 0  0 . 7 5 0 0 −0.5260 −0.4381 −0.3590];
ychord ( 1 , 2 8 : 3 6 ) = [  0 . 0 2 1 2  0 . 0 1 5 5  0 . 0 1 0 5  0 . 0 0 5 6  0 . 0 0 0 6 −0.0056 −0.0106 −0.0155 −0.0205];
ychord ( 2 , 2 8 : 3 6 ) = [ −0.2877 −0.2102 −0.1497 −0.0927 −0.0364  0 . 0 3 1 4  0 . 0 8 7 7  0 . 1 4 4 7  0 . 2 0 5 2 ] ;
ychord ( 1 , 3 7 : 4 5 ) = [ −0.0262 −0.0305 −0.0348 −0.0391 −0.4234 −0.4277 −0.4320 −0.4363 −0.4420];
ychord ( 2 , 3 7 : 4 5 ) = [ 0 . 2 8 2 7 0 . 3 5 4 0 0 . 4 3 3 1 0 . 5 2 1 0 −0.7500 −0.7500 −0.7500 −0.7452 −0.6677];
ychord ( 1 , 4 6 : 5 4 ) = [ −0.4470 −0.4519 −0.4569 −0.4631 −0.4681 −0.4730 −0.4780 −0.4837 −0.4880];
ychord ( 2 , 4 6 : 5 4 ) = [ −0.6072 −0.5502 −0.4939 −0.4261 −0.3698 −0.3128 −0.2523 −0.1748 −0.1035];
ychord ( 1 , 5 5 : 6 3 ) = [ −0.4923 −0.4966 −0.9057 −0.9063 −0.9069 −0.9075 −0.9075 −0.9069 −0.9063];
ychord ( 2 , 5 5 : 6 3 ) = [ −0.0244 0 . 0 6 3 5 0 . 4 9 5 1 0 . 7 5 0 0 0 . 7 5 0 0 0 . 7 5 0 0 0 . 7 5 0 0 0 . 7 1 0 0 0 . 6 2 0 0 ] ;
ychord ( 1 , 6 4 : 6 4 ) = [ −0.9057];
ychord ( 2 , 6 4 : 6 4 ) = [  0 . 4 8 0 0 ] ;
save fabchord . mat xchord ychord ;

disp ( 'File fabchord.mat has been created successfully.' ) ;

disp ( ' ' ) ;

B.2.5 fabsetsolidangle.m

function [OMEGA, r h o g r i d , z e t g r i d ] = f a b s e t s o l i d a n g l e
% f u n c t i o n [OMEGA, r h o g r i d , z e t g r i d ] = f a b s e t s o l i d a n g l e

%

% ’ f a b s e t s o l i d a n g l e .m’ c a l c u l a t e s the i n d i v i d u a l bolometers s o l i d angles and s t o r e s
% them to f i l e s c a l l e d ’ f a b s o l i d a n g l e #.mat ’ in the current working d i r e c t o r y .

%

% The ’ f a b s o l i d a n g l e #.mat ’ f i l e s are used by the r o u t i n e ’ f a b s o l i d a n g l e .m’ which helps
% to c a l c u l a t e the T−matrix ( i n i t i a t e d by ’ f a b s e t t .m’ )

%

% This program c a l l s ’ fab3d .m’

%

% History :

%

% Matthias Anton , CRPP EPFL, 2 9 of May, 1 9 9 5 ( f i r s t r e l e a s e )
% Jan Mlynar , CRPP EPFL, 2 0 of August , 1 9 9 8 ( adaption f o r bolometry )
% C h r i s t i a n S c h l a t t e r , TPIV CRPP EPFL , February 2 0 0 2 ( comments and renaming )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
global Kb1 Kb2 Kb3 Kb4 Kd1 Kd2 Kd3 Kd4
debut=cputime ;
disp ( ' ' )

disp ( 'Solid angle calculator for bolometry.' ) ;
disp ( 'Final release 1.4' )
disp ( 'Latest modification on 2002/02/07.' )

disp ( ' ' ) ;
% Detector parameters
cw=1;  % Detector numbers cw = 1 : c l o c k w i s e cw = 0 : ccw
fans = [ 1 1 1 1 1 1 1 1 ] ;  % Camera switch
vangle = [ 9 0 0 0 0 0 0 0 − 9 0 ] ; % Angle of d e t e c t o r s u r f a c e normal
xpos = [ 8 8 1 2 3 . 5 1 2 3 . 5 1 2 3 . 5 1 2 3 . 5 1 2 3 . 5 1 2 3 . 5 8 8 ] ;
% x p o s i t i o n of the diaphragmas in [ cm ]
ypos =[81.5 45.5 45.5 −.25 −.25 −46 −46 −81.5];
% y p o s i t i o n of the diaphragmas in [ cm ]
% P o s i t i o n of the d e t e c t o r s
rdet = [ 0 . 8 6 2 , 0 . 8 6 7 , 0 . 8 7 2 1 , 0 . 8 7 7 1 , 0 . 8 8 2 9 , 0 . 8 8 7 9 , 0 . 8 9 3 , 0 . 8 9 8 , …

1 . 2 7 8 , 1 . 2 8 0 6 , 1 . 2 8 3 3 , 1 . 2 8 6 1 , 1 . 2 8 8 3 , 1 . 2 8 9 3 , 1 . 2 9 0 2 , 1 . 2 9 1 2 ] ;

rdet (17:24)= f l i p l r ( rdet ( 9 : 1 6 ) ) ;
rdet (25:40)= rdet ( 9 : 2 4 ) ;
rdet (41:56)= rdet ( 9 : 2 4 ) ;
rdet (57:64)= f l i p l r ( rdet ( 1 : 8 ) ) ;
zdet = [ 0 . 9 0 5 7 , 0 . 9 0 6 3 , 0 . 9 0 6 9 , 0 . 9 0 7 5 ] ;
zdet (5:8)= f l i p l r ( zdet ( 1 : 4 ) ) ;
zdet ( 9 : 2 4 ) = [ 0 . 4 9 1 6 , 0 . 4 8 7 3 , 0 . 4 8 3 , 0 . 4 7 8 7 , 0 . 4 7 3 , 0 . 4 6 8 , 0 . 4 6 3 1 , 0 . 4 5 8 1 , …

0 . 4 5 1 9 , 0 . 4 4 6 9 , 0 . 4 4 2 , 0 . 4 3 7 , 0 . 4 3 1 3 , 0 . 4 2 7 , 0 . 4 2 2 7 , 0 . 4 1 8 4 ] ;

zdet (25:40)= zdet ( 9 : 2 4 ) − 0 . 4 5 7 5 ;
zdet (41:56)= zdet ( 9 : 2 4 ) − 0 . 9 1 5 ;
zdet (57:64)=− zdet ( 1 : 8 ) ;
xdet=rdet ∗ 1 0 0 ; ydet=zdet ∗100;
d1 =1.5/10;  % Detector width  in cm
d2=4/10;  % Detector length  in cm
b1 = [ 2 . 6 2 . 2 2 . 2 2 . 2 2 . 2 2 . 2 2 . 2 2 . 6 ] / 5 ;
% Aperture width  in cm ( pol . )
b2 = [ 1 0 8 8 8 8 8 8 1 0 ] / 5 ;
% Aperture length  in cm ( tor . )
% C a l c u l a t i o n of aperture and d e t e c t o r p o s i t i o n s
nact=sum( fans ) ;
i a c t=find ( fans ) ;
ndet =8;
ncam=8;

% Aperture 1 ( p o l o i d a l l y l i m i t i n g )
xap=ones ( ndet , 1 ) ∗ xpos ;
xap=xap ( : ) ’ ;
yap=ones ( ndet , 1 ) ∗ ypos ;
yap=yap ( : ) ’ ;
% Detectors
dxd1=xdet−xap ;
dyd1=ydet−yap ;
i v e r t = [ 1 : 8 , 5 7 : 6 4 ] ’ ;
i h o r i = [ 9 : 5 6 ] ’ ;
dxd2=dxd1 ; dyd2=dyd1 ; % Supp . da = da2 i . e . diaphragma−array d i s t a n c e i s
% equal in p o l o i d a l ( da ) and t o r o i d a l ( da2 ) d i r e c t i o n s
xap2=xdet−dxd2 ;
yap2=ydet−dyd2 ;
xyw=ones ( ndet , 1 ) ∗ b1 /2;
xyw=xyw ( : ) ’ ;
idxw ( i v e r t )=ones ( size ( i v e r t ) ) ;
idyw ( i h o r i )=ones ( size ( i h o r i ) ) ;
idyw (64)=0;
dxw=xyw ;
dyw=xyw ;
dzw=ones ( ndet , 1 ) ∗ b2 ( i a c t ) / 2 ;
dzw=dzw ( : ) ’ ;
% x , y : p o l o i d a l c r o s s s e c t i o n ; z : ’ t o r o i d a l ’ coordinate
% Transform to the ’ format ’ used by ’ fab3d .m’ :
% 1 , 2 : midpoints on the s t r a i g h t l i n e s which l i m i t the d e t e c t o r p o l o i d a l l y
% 3 , 4 ; midpoints on the s t r a i g h t l i n e s which l i m i t the d e t e c t o r t o r o i d a l l y
Kb1=[xap−dxw ; yap−dyw ; zeros ( size ( xap ) ) ] ;
Kb2=[xap+dxw ; yap+dyw ; zeros ( size ( xap ) ) ] ;
Kb3=[xap2 ; yap2;−dzw ] ;
Kb4=[xap2 ; yap2 ; dzw ] ;
Kd1=[ xdet−idxw∗d1 /2; ydet−idyw∗d1 /2; zeros ( size ( xap ) ) ] ;
Kd2=[ xdet+idxw∗d1 /2; ydet+idyw∗d1 /2; zeros ( size ( xap ) ) ] ;
Kd3=[ xdet ; ydet ;−d2∗ ones ( size ( xdet ) ) / 2 ] ;
Kd4=[ xdet ; ydet ; d2∗ ones ( size ( xdet ) ) / 2 ] ;
% Get the d e t e c t o r names

APPENDIX B. THE PROGRAM CODE 66

detn = [ ] ;
for j =1:length ( fans )
i f fans ( j )
detn =[ detn , ( j −1)∗8+1: j ∗ 8 ] ;

end end

% Calculate the s o l i d angle f o r the whole g r i d and every d e t e c t o r
for kk=1:length ( xap )
disp ( '' ) ;

    disp ( [ '  Calculation of the solid angle for detector ' , int2str ( kk ) ] ) ;

[ omega , r h o g r i d , z e t g r i d , dV ] = …
fab3d (Kb1 ( : , kk ) , Kb2 ( : , kk ) , Kb3 ( : , kk ) , Kb4 ( : , kk ) , …
Kd1 ( : , kk ) , Kd2 ( : , kk ) , Kd3 ( : , kk ) , Kd4 ( : , kk ) , idxw ( kk ) ) ;
i c i=find ( omega ) ;
rhodum=ones ( size ( z e t g r i d ’ ) ) ∗ r h o g r i d ;
zetdum=z e t g r i d ’ ∗ ones ( size ( r h o g r i d ) ) ;
r h o i c i=rhodum ( i c i ) ;
z e t i c i=zetdum ( i c i ) ;
omeici=omega ( i c i ) ;
i d e t e c =1;
eval ( [ 'save fabsolidangle_' , int2str ( detn ( kk ) ) , . . . . . .

         '.mat ici rho_grid zet_grid rhoici zetici omeici dV i_detec' ] ) ;

end

OMEGA= [ ] ;
disp ( ' ' ) ;

disp ( 'Successfully finished' ) ;

disp ( ' ' ) ;
time end = ( cputime − debut ) / 6 0 ;
disp ( [ 'CPU runtime: ' num2str( time end ) , ' min.' ] ) ;
disp ( ' ' ) ;

B.2.6 fab3d.m

function [ omega , r h o g r i d , z e t g r i d , dV ] = …

fab3d (Kb1 , Kb2 , Kb3 , Kb4 , Kd1 , Kd2 , Kd3 , Kd4 , i v e r t )
% f u n c t i o n [ omega , r h o g r i d , z e t g r i d , dV ] = …
%  fab3d (Kb1 , Kb2 , Kb3 , Kb4 , Kd1 , Kd2 , Kd3 , Kd4 , i v e r t )

%

% C a l c u l a t e s a 2D−matrix omega from a 3D g r i d definded i n s i d e
% ’ f a b s e t a n g l e .m’ .
% The g r i d f i l l s approximately a 4 0 cm t h i c k p o l o i d a l s l i c e
% of the TCV vacuum v e s s e l ( ’ t h i c k ’ : in t o r o i d a l d i r e c t i o n )

%

% input data : Kb1 . . 4 :  midpoints of the edges of the aperture
% Kd1 . . 4 : midpoints of the edges of the d e t e c t o r
%  each K . . has three components :
%  K. . ( 1 ) : r a d i a l rho
%  K. . ( 2 ) : v e r t i c a l  zet
%  K. . ( 3 ) : t o r o i d a l  tee
%  K. . 1 & 2 : midpoints of edge ’ l i n e s ’
%  in rho−zet−plane
%  K. . 3 & 4 : midpoints of edges
%  in tee−zet−plane
% i v e r t : a f l a g , determines i f the d e t e c t o r ’ l o o k s ’
%  h o r i z o n t a l l y or v e r t i c a l l y
% output data : omega :  The d e t e c t o r s s o l i d angles
%  r h o g r i d : The r a d i a l c o o r d i n a t e s .
%  z e t g r i d : The v e r t i c a l c o o r d i n a t e s .
% dV : a p i x e l s volume .

%

% This program c a l l s ’ f a b p r o j e c t .m’

%

% History :

%

% Matthias Anton , CRPP EPFL, 2 9 of May, 1 9 9 5 ( f i r s t r e l e a s e )
% C h r i s t i a n S c h l a t t e r , TPIV CRPP EPFL , February 2 0 0 2 ( comments and renaming )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )

% parameter : d e f i n e the 3D g r i d
% rho : r a d i a l
% zet : v e r t i c a l
% tee : t o r o i d a l coordinate

Rho=88;
brho =2∗30;
bzet =2∗90;
btee =2∗20;
drho =0.5;

dzet =0.5;
dtee =0.5;
dV=dtee ∗ dzet ∗drho ;
r h o g r i d=Rho−(brho−drho ) / 2 : drho : Rho+(brho−drho ) / 2 ;
z e t g r i d=−(bzet−dzet ) / 2 : dzet : ( bzet−dzet ) / 2 ;
t e e g r i d =0: dtee : ( btee−dtee ) / 2 ;
% get four points of aperture ( midpoint of each s i d e )
% 1 , 2 : in the rho−zet plane
% 3 , 4 : in the tee−zet plane
rhob1=Kb1 ( 1 ) ;
zetb1=Kb1 ( 2 ) ;
teeb1=Kb1 ( 3 ) ;
rhob2=Kb2 ( 1 ) ;
zetb2=Kb2 ( 2 ) ;
teeb2=Kb2 ( 3 ) ;
rhob3=Kb3 ( 1 ) ;
zetb3=Kb3 ( 2 ) ;
teeb3=Kb3 ( 3 ) ;
rhob4=Kb4 ( 1 ) ;
zetb4=Kb4 ( 2 ) ;
teeb4=Kb4 ( 3 ) ;
% get four points of d e t e c t o r ( midpoint of each s i d e )
% 1 , 2 : in the rho−zet plane
% 3 , 4 : in the tee−zet plane
rhod1=Kd1 ( 1 ) ;
zetd1=Kd1 ( 2 ) ;
teed1=Kd1 ( 3 ) ;
rhod2=Kd2 ( 1 ) ;
zetd2=Kd2 ( 2 ) ;
teed2=Kd2 ( 3 ) ;
rhod3=Kd3 ( 1 ) ;
zetd3=Kd3 ( 2 ) ;
teed3=Kd3 ( 3 ) ;
rhod4=Kd4 ( 1 ) ;
zetd4=Kd4 ( 2 ) ;
teed4=Kd4 ( 3 ) ;

% i n i t omega
omega=zeros ( length ( z e t g r i d ) , length ( r h o g r i d ) ) ;
% loop over t o r o i d a l coordinate
for i t =1:length ( t e e g r i d ) ;
tee=t e e g r i d ( i t ) ;
r h o s h i f t=Rho−sqrt (Rhoˆ 2−tee ˆ 2 ) ;
rho=ones ( size ( z e t g r i d ’ ) ) ∗ ( r h o g r i d−r h o s h i f t ) ;
zet=z e t g r i d ’ ∗ ones ( size ( r h o g r i d ) ) ;
tee=tee ∗ ones ( size ( zet ) ) ;
i f i v e r t
[ rhol1 , z e t l 1 , rhol2 , z e t l 2 ] = …
f a b p r o j e c t ( rhob1 , zetb1 , rhob2 , zetb2 , rhod1 , zetd1 , rhod2 , zetd2 , rho , zet ) ;
[ t e e l 3 , z e t l 3 , t e e l 4 , z e t l 4 ] = …
f a b p r o j e c t ( teeb3 , zetb3 , teeb4 , zetb4 , teed3 , zetd3 , teed4 , zetd4 , tee , zet ) ;
die=find ( rhol1 ˜ =99999 & t e e l 3 ˜ =99999);
d i s r h o =( rhol1 ( die )+ rhol2 ( die ))/2− rho ( die ) ;
d i s z e t =( z e t l 1 ( die )+ z e t l 2 ( die )+ z e t l 3 ( die )+ z e t l 4 ( die ))/4− zet ( die ) ;
d i s t e e =( t e e l 3 ( die )+ t e e l 4 ( die ))/2− tee ( die ) ;
s u r f r h o =( z e t l 1 ( die )− z e t l 2 ( die ) ) . ∗ ( t e e l 3 ( die )− t e e l 4 ( die ) ) ;
s u r f z e t=−(rhol1 ( die )− rhol2 ( die ) ) . ∗ ( t e e l 3 ( die )− t e e l 4 ( die ) ) ;
s u r f t e e =( rhol1 ( die )− rhol2 ( die ) ) . ∗ ( z e t l 3 ( die )− z e t l 4 ( die ) ) ;

else

[ z e t l 1 , rhol1 , z e t l 2 , rhol2 ] = ….
f a b p r o j e c t ( zetb1 , rhob1 , zetb2 , rhob2 , zetd1 , rhod1 , zetd2 , rhod2 , zet , rho ) ;
[ t e e l 3 , rhol3 , t e e l 4 , rhol4 ] = ….
f a b p r o j e c t ( teeb3 , rhob3 , teeb4 , rhob4 , teed3 , rhod3 , teed4 , rhod4 , tee , rho ) ;
die=find ( z e t l 1 ˜ =99999 & t e e l 3 ˜ =99999);
d i s r h o =( rhol1 ( die )+ rhol2 ( die )+ rhol3 ( die )+ rhol4 ( die ))/4− rho ( die ) ;
d i s z e t =( z e t l 1 ( die )+ z e t l 2 ( die ))/2− zet ( die ) ;
d i s t e e =( t e e l 3 ( die )+ t e e l 4 ( die ))/2− tee ( die ) ;
s u r f r h o =( z e t l 1 ( die )− z e t l 2 ( die ) ) . ∗ ( t e e l 3 ( die )− t e e l 4 ( die ) ) ;
s u r f z e t=−(rhol1 ( die )− rhol2 ( die ) ) . ∗ ( t e e l 3 ( die )− t e e l 4 ( die ) ) ;

APPENDIX B. THE PROGRAM CODE 70

s u r f t e e=−(z e t l 1 ( die )− z e t l 2 ( die ) ) . ∗ ( rhol3 ( die )− rhol4 ( die ) ) ;

end

dd=( d i s r h o .ˆ 2+ d i s z e t .ˆ 2+ d i s t e e . ˆ 2 ) . ˆ 0 . 5 ;
s u r f n=abs ( s u r f r h o . ∗ d i s r h o + …
s u r f z e t . ∗ d i s z e t + …
s u r f t e e . ∗ d i s t e e ) . / dd ;
omega ( die )=omega ( die ) …
+abs ( s u r f r h o . ∗ d i s r h o + …
s u r f z e t . ∗ d i s z e t + …
s u r f t e e . ∗ d i s t e e ) . / dd . ˆ 3 ;
i f i t ==1;
omega0=omega ;

end end

omega=2∗omega−omega0 ;

B.2.7 fabproject.m

function [ xl1 , yl1 , xl2 , yl2 ] = …
f a b p r o j e c t ( xb1 , yb1 , xb2 , yb2 , xd1 , yd1 , xd2 , yd2 , xi , yi )
% f u n c t i o n [ xl1 , yl1 , xl2 , yl2 ] = …
%  f a b p r o j e c t ( xb1 , yb1 , xb2 , yb2 , xd1 , yd1 , xd2 , yd2 , xi , yi )

%

% C a l c u l a t e s the p r o j e c t i o n of two points xb1 , xb2 , yb1 , yb2 on a l i n e
% defined by xd1 , yd1 , xd2 , yd2 . The point of p r o j e c t i o n i s xi , yi .

%

% s i z e s : xb1 , yb1 , xb2 , yb2 , xd1 , yd1 , xd2 , yd2  1x1
%  xi , yi  a r b i t r a r y
%  xl1 , yl1 , xl2 , yl2  same s i z e as xi , yi

%

% This program i s a subroutine of ’ fab3d .m’

%

% History :

%

% Matthias Anton , CRPP EPFL, 2 9 of May, 1 9 9 5 ( f i r s t r e l e a s e )
% C h r i s t i a n S c h l a t t e r , TPIV CRPP EPFL , February 2 0 0 2 ( comments and renaming )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
xl1=zeros ( size ( xi ) ) ;

xl2=zeros ( size ( xi ) ) ;
yl1=zeros ( size ( xi ) ) ;
yl2=zeros ( size ( xi ) ) ;
p1=( yi ∗( xb1−xd1)− xi ∗( yb1−yd1)+xd1∗yb1−yd1∗xb1 ) . / ….
( ( yb1−yi )∗( xd1−xd2)−(xb1−xi )∗( yd1−yd2 ) ) ;
p2=( yi ∗( xb2−xd1)− xi ∗( yb2−yd1)+xd1∗yb2−yd1∗xb2 ) . / ….
( ( yb2−yi )∗( xd1−xd2)−(xb2−xi )∗( yd1−yd2 ) ) ;
xb1d=xd1+p1 ∗( xd2−xd1 ) ;
yb1d=yd1+p1 ∗( yd2−yd1 ) ;
xb2d=xd1+p2 ∗( xd2−xd1 ) ;
yb2d=yd1+p2 ∗( yd2−yd1 ) ;
i i a = find ( ( p1<0 & p2 < 0 ) | ( p1>1 & p2 > 1 ) ) ;
i i b = find ( p1<=1 & p1>=0 & p2>=0 & p2<=1 );
i i c = find ( ( p1<0 & p2 > 1 ) | ( p2<0 & p1 > 1 ) ) ;
i i d = find ( p1>=0 & p1<=1 & p2 >1 );
i i e = find ( p1>=0 & p1<=1 & p2 <0 );
i i f = find ( p2>=0 & p2<=1 & p1 >1 );
i i g = find ( p2>=0 & p2<=1 & p1 <0 );
i f ˜ isempty ( i i a ) ;
xl1 ( i i a )= xi ( i i a ) ;
xl2 ( i i a )= xi ( i i a ) ;
yl1 ( i i a )= yi ( i i a ) ;
yl2 ( i i a )= yi ( i i a ) ;

end

i f ˜ isempty ( i i a ) ;
xl1 ( i i a )=99999∗ ones ( size ( i i a ) ) ;
xl2 ( i i a )=99999∗ ones ( size ( i i a ) ) ;
yl1 ( i i a )=99999∗ ones ( size ( i i a ) ) ;
yl2 ( i i a )=99999∗ ones ( size ( i i a ) ) ;

end

i f ˜ isempty ( i i b ) ;
xl1 ( i i b )=xb1d ( i i b ) ;
xl2 ( i i b )=xb2d ( i i b ) ;
yl1 ( i i b )=yb1d ( i i b ) ;
yl2 ( i i b )=yb2d ( i i b ) ;

end

i f ˜ isempty ( i i c ) ;
xl1 ( i i c )=xd1∗ ones ( size ( i i c ) ) ;
xl2 ( i i c )=xd2∗ ones ( size ( i i c ) ) ;
yl1 ( i i c )=yd1∗ ones ( size ( i i c ) ) ;
yl2 ( i i c )=yd2∗ ones ( size ( i i c ) ) ;

end

i f ˜ isempty ( i i d ) ;
xl1 ( i i d )=xd2∗ ones ( size ( i i d ) ) ;
xl2 ( i i d )=xb1d ( i i d ) ;
yl1 ( i i d )=yd2∗ ones ( size ( i i d ) ) ;
yl2 ( i i d )=yb1d ( i i d ) ;

end

i f ˜ isempty ( i i e ) ;
xl1 ( i i e )=xd1∗ ones ( size ( i i e ) ) ;
xl2 ( i i e )=xb1d ( i i e ) ;
yl1 ( i i e )=yd1∗ ones ( size ( i i e ) ) ;
yl2 ( i i e )=yb1d ( i i e ) ;

end

i f ˜ isempty ( i i f ) ;
xl1 ( i i f )=xd2∗ ones ( size ( i i f ) ) ;
xl2 ( i i f )=xb2d ( i i f ) ;
yl1 ( i i f )=yd2∗ ones ( size ( i i f ) ) ;
yl2 ( i i f )=yb2d ( i i f ) ;

end

i f ˜ isempty ( i i g ) ;
xl1 ( i i g )=xd1∗ ones ( size ( i i g ) ) ;
xl2 ( i i g )=xb2d ( i i g ) ;
yl1 ( i i g )=yd1∗ ones ( size ( i i g ) ) ;
yl2 ( i i g )=yb2d ( i i g ) ;

end

B.2.8 fabsett.m

function f a b s e t t

% f u n c t i o n f a b s e t t

%

% C a l c u l a t i o n of the T matrix f o r a r e c t a n g u l a r g r i d using p r e c a l c u l a t e d
% matrices of s o l i d angles f o r a l l d e t e c t o r s and a g r i d of 0 . 5 x0 . 5 x0 . 5cmˆ 3
% The s o l i d angle matrices has to be c a l l e d ’ f a b a n g l e #.mat ’ and need to be
% stored in a s u b d i r e c t o r y c a l l e d \ f a b s o l i d a n g l e \.

%

% This uses the 3 dimensional omgrid algorithm ( fabangle .m)
% This i s the v e r s i o n f o r bolometry .

%

% The v a r i a b l e s nx ny xmin xmax ymin ymax dx dy xmesh ymesh numdet T
% are stored to the f i l e fabt . mat f o r use with fabcat .m.

%

% History :

%

% C h r i s t i a n S c h l a t t e r , TPIV CRPP EPFL , November 2001

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
disp ( ' ' )

disp ( 'FABCAT T-matrix calculator for bolometry.' ) ;
disp ( 'Final release 1.4' )
disp ( 'Latest modification on 2002/02/07.' )

disp ( ' ' ) ;
% Tokamak setup parameters f o r p l a c i n g the r e c t a n g u l a r g r i d .
% Values correspond to the d e f a u l t values in the btomo package .
% All values given in [ cm ]
r 0 = 8 8 ;  % r a d i a l p o s i t i o n of the r e c t a n g u l a r g r i d c e n t e r
z 0 = 0 ;  % z−coordinate of the r e c t a n g u l a r g r i d c e n t e r
w r = 5 5 ;  % h o r i z o n t a l width of the g r i d
w z = 1 5 4 ;  % v e r t i c a l width of the g r i d
nx = 1 0 ;  % h o i z o n t a l number of p i x e l s
% Camera setup ( bolometry d e t e c t o r s )
f cam = [ 1 1 1 1 1 1 1 1 ] ;
xmin=r 0−w r /2;
xmax=r 0+w r /2;
ymin=z 0−w z /2;
ymax=z 0+w z /2;
dx=(xmax−xmin )/ nx ;
dy=dx ;
ny = c e i l ( ( ymax−ymin )/ dy ) ;

APPENDIX B. THE PROGRAM CODE 74

ymin = z 0−ny/2∗dy ;
ymax = z 0+ny/2∗dy ;
xmesh = xmin+dx /2: dx : xmax−dx /2;
ymesh = ymin+dy /2: dy : ymax−dy /2;
[T, numdet ] = fabangle ( f cam , xmin , xmax , ymin , ymax , nx , ny ) ;
save fabt . mat nx ny xmin xmax ymin ymax dx dy xmesh ymesh numdet T;

disp ( 'File fabt.mat successfully created.' )

disp ( ' ' ) ;

B.2.9 fabangle.m

function [ T , numdet ] = fabangle ( fans , xmin , xmax , ymin , ymax , nx , ny )
% f u n c t i o n [ T , numdet ] = fabangle ( fans , xmin , xmax , ymin , ymax , nx , ny )

%

% C a l c u l a t i o n of the T matrix f o r a r e c t a n g u l a r g r i d using p r e c a l c u l a t e d
% matrices of s o l i d angles f o r a l l d e t e c t o r s and a g r i d of 0 . 5 x0 . 5 x0 . 5cmˆ 3
% matrices are stored in f a b a n g l e #.mat f i l e s l o c a t e d in the s u b d i r e c t o r y
% \ f a b s o l i d a n g l e \.

%

% History :

%

% Matthias Anton , CRPP EPFL, 3 0 of May, 1 9 9 5 ( f i r s t r e l e a s e )
% Matthias Anton , CRPP EPFL, 2 3 of November , 1 9 9 5
% Iwan J e r j e n , CRPP EPFL, 1 2 of July , 1 9 9 8
% Jan Mlynar , CRPP EPFL, 2 0 of August , 1 9 9 8 ( adaption f o r bolometry )

%

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
debut=cputime ;
detn = [ ] ;
for j =1:length ( fans )
i f fans ( j )
detn =[ detn , ( j −1)∗8+1: j ∗ 8 ] ;

end end

ndet=length ( detn ) ;
dx=(xmax−xmin )/ nx ;

APPENDIX B. THE PROGRAM CODE 75

dy=(ymax−ymin )/ ny ;
xgrid=xmin+dx /2: dx : xmax−dx /2;
ygrid=ymin+dy /2: dy : ymax−dy /2;
xpix=ones ( size ( ygrid ’ ) ) ∗ xgrid ;
ypix=ygrid ’ ∗ ones ( size ( xgrid ) ) ;
xpix=xpix ( : ) ;
ypix=ypix ( : ) ;
T=zeros ( ndet , nx∗ny ) ;
for k=1: ndet

  eval ( [ 'load fabsolidangle/fabsolidangle_' . . .

, int2str ( detn ( k ) ) ] )
for l = 1 : nx∗ny
drin=find ( r h o i c i>=xpix ( l )−dx/2 & r h o i c i<=xpix ( l )+dx / 2 ….
& z e t i c i>=ypix ( l )−dy/2 & z e t i c i<=ypix ( l )+dy / 2 ) ;
T( k , l )=sum( omeici ( drin ))∗dV;

end end

load fabchord
[ Th , numh]= fabstandard (100.∗ xchord , 1 0 0 . ∗ ychord , xmin , xmax , ymin , ymax , nx , ny ) ;
for i =1:64
l i n s t=find (numh=i ) ;
i f ˜ isempty ( l i n s t )
leng=sum(Th( l i n s t , : ) ) ;

else

leng =0;

end

viv=sum(T( i , : ) ) ;
i f viv˜ =0
cor=leng / viv ;

APPENDIX B. THE PROGRAM CODE 76

else

cor =0;

end

T( i , : ) =T( i , : ) ∗ cor ;

end

numdet=find (sum(T ’ ) ) ;
T=T( numdet , : ) ;
numdet=detn ( numdet ) ;
time end = cputime − debut ;
disp ( [ 'CPU runtime: ' num2str( time end ) , ' s.' ] ) ;
disp ( '  ' ) ;

B.2.10 fabstandard.m

function [TT, numdet ] = fabstandard ( xchord , ychord , …
xmin , xmax , ymin , ymax , nx , ny ) ;
% f u n c t i o n [TT, numdet ] = fabstandard ( xchord , ychord , …
%  xmin , xmax , ymin , ymax , nx , ny ) ;

%

% A f a s t algorithm to c a l c u l a t e the l e n g t h s of the chords given by
% xchord , ychord in p i x e l s of a g r i d s p e c i f i e d by the other inputs .

%

% Determination only based on geometry

%

% Inputs

%

%  xchord , ychord : endpoints of l i n e s of s i g h t , s i z e [ 2 x nl ]
%  xmin … ymax : c o r n e r s of p i x e l g r i d
% nx , ny : number of p i x e l s h o r i z o n t a l , v e r t i c a l

%

% Output

%

% TT: t r a n s f e r m a t r i x [ length ( numdet ) x nx∗ny ] ,
%  TT( l , i ) i s the length of chord l in p i x e l i
%  numdet : numbers of ’ a c t i v e ’ l i n e s of s i g h t , u s u a l l y
%  length ( numdet ) < nl

%

% History :

%

% Matthias Anton , CRPP EPFL , 9 of August 1 9 9 4 ( f i r s t r e l e a s e )
% Matthias Anton , CRPP EPFL, 2 nd of December 1994

APPENDIX B. THE PROGRAM CODE 77 %

% Final r e l e a s e v e r s i o n 1 . 4 ( 2 0 0 2 / 0 2 / 0 7 )
[dummy, nl ]= size ( xchord ) ;
dx=(xmax−xmin )/ nx ;
dy=(ymax−ymin )/ ny ;
xgrid=xmin : dx : xmax ;
ygrid=ymin : dy : ymax ;
xpix=xmin+dx /2: dx : xmax−dx /2;
ypix=ymin+dy /2: dy : xmax−dy /2;
numpix=reshape ( 1 : nx∗ny , ny , nx ) ;
for k=1: nl
c=polyfit ( xchord ( : , k ) , ychord ( : , k ) , 1 ) ;
m( k)=c ( 1 ) ; b( k)=c ( 2 ) ;

end

% c r o s s i n g with v e r t i c a l l i n e s of g r i d
Ysec=zeros ( nl , nx+1);
X=ones ( nl , 1 ) ∗ xgrid ;
Ysec=m’ ∗ xgrid+b ’ ∗ ones ( 1 , nx+1);
% c r o s s i n g s with h o r i z o n t a l l i n e s of g r i d
Xsec=zeros ( nl , ny+1);
Y=ones ( nl , 1 ) ∗ ygrid ;
Xsec=(ones ( nl , 1 ) ∗ ygrid−b ’ ∗ ones ( 1 , ny +1))./(m’ ∗ ones ( 1 , ny +1));
% matrices with x and y c o o r d i n a t e s of a l l c r o s s i n g s with meshgrid
XX=[X, Xsec ] ;
YY=[Ysec ,Y ] ;
% s o r t i n g in ascending order of X
for k=1: nl
[XX( k , : ) , ind ]= sort (XX( k , : ) ) ;
YY( k , : ) =YY( k , ind ) ;

end

dX=[ d i f f (XX’ ) ’ zeros ( nl , 1 ) ] ;
dY=[ d i f f (YY’ ) ’ zeros ( nl , 1 ) ] ;

APPENDIX B. THE PROGRAM CODE 78

% t r y i n g to f i n d out the p i x e l numbers
XXnum=zeros ( size (XX) ) ;
YYnum=zeros ( size (YY) ) ;
for k=1:nx
hx=find (XX+dX/2>xgrid ( k ) & XX+dX/2<xgrid ( k+1) & dX>0);
XXnum( hx)=k∗ ones ( size ( hx ) ) ;

end

for k=1:ny
hy=find ( (YY+dY/2>ygrid ( k ) & YY + dY/2 < ygrid ( k+1)) & dX>0 );
YYnum( hy)=k∗ ones ( size ( hy ) ) ;

end

s e g l e n g t h=sqrt (dX.ˆ 2+dY . ˆ 2 ) ;
TT=zeros ( nl , nx∗ny ) ;
for k=1: nl
c o l s=find ( (YYnum( k ,:)˜ =0)&(XXnum( k , : ) ˜ = 0 ) & …
( s e g l e n g t h ( k,:)<= sqrt ( dxˆ 2+dy ˆ 2 ) ) ) ;
dummy=diag ( numpix (YYnum( k , c o l s ) ,XXnum( k , c o l s ) ) ) ;
i f ˜ isempty (dummy)
TT( k ,dummy(:))= s e g l e n g t h ( k , c o l s ) ;

end end

numdet=find (sum(TT’ ) ) ;
TT=TT( numdet , : ) ;

B.3 Tools

B.3.1 fabmov.m

function bolmov ( shot , wperclev )
% bolmov ( shot , l e v e l ) / without argument runs l a s t shot , l e v e l =100

%

% Loads r e s u l t s of f a s t SXR tomography . Uses BLACK AND WHITE contours f o r
% d i s p l a y i n g emission p r o f i l e s . Shows plasma p o s i t i o n , allows MOVIE run .

% By d e f a u l t the emission l e v e l i s 2 0 kW.mˆ (−3) but t h i s can be changed
% by the second f u n c t i o n parameter .
% Jan Mlynar Sept 2001
% v e r s i o n bolometry … January 2002
global bm maincat bm movieaxis bm shotnum bm comment bm framenum bm chisq …
bm lambda bm Ptot bm Pup bm Pbot bm play bm repeat bm speed …
bm delchans bm delchans2
global shotnum t s t e p s comment g bolo Ptot lambda rpix zpix tvec c h i s q …
rmag zmag moma maxlev c l e v s wperc Pup Pbot delchans
i f nargin<1
shot =21518;

end

i f isnumeric ( shot )
a c t i o n='start' ;
shotnum=shot ;
i f nargin<2
wperc =20000; % W.mˆ (−3) per contour

else

wperc=wperclev ;

end else

a c t i o n=shot ;

end

switch a c t i o n

case 'start'

bm maincat=bolmovgui ;
zoom on ;
bm movieaxis=f i n d o b j ( bm maincat , 'Tag' , 'movieaxes' ) ;
bm shotnum=f i n d o b j ( bm maincat , 'Tag' , 'shotnum' ) ;
bm comment=f i n d o b j ( bm maincat , 'Tag' , 'comment' ) ;
bm framenum=f i n d o b j ( bm maincat , 'Tag' , 'framenum' ) ;
bm chisq=f i n d o b j ( bm maincat , 'Tag' , 'chisq' ) ;
bm lambda=f i n d o b j ( bm maincat , 'Tag' , 'lambda' ) ;
bm Ptot=f i n d o b j ( bm maincat , 'Tag' , 'Ptot' ) ;
bm Pup=f i n d o b j ( bm maincat , 'Tag' , 'Pup' ) ;
bm Pbot=f i n d o b j ( bm maincat , 'Tag' , 'Pbot' ) ;
bm play=f i n d o b j ( bm maincat , 'Tag' , 'play' ) ;
bm repeat=f i n d o b j ( bm maincat , 'Tag' , 'repeat' ) ;
bm speed=f i n d o b j ( bm maincat , 'Tag' , 'speed' ) ;
bm delchans=f i n d o b j ( bm maincat , 'Tag' , 'delchans' ) ;

bm delchans2=f i n d o b j ( bm maincat , 'Tag' , 'delchans2' ) ;

axes ( bm movieaxis ) ;
axis equal ;
axis ( [ . 6 1 . 2 −. 8 . 8 ] ) ;

set ( bm shotnum , 'String' , int2str ( shotnum ) ) ;

bolmov ( 'shotnum' ) ;

return

case 'shotnum'
set (bm comment , 'String' , '...wait please, getting inversions...' ) ;

cla drawnow

shotnum=eval ( get ( bm shotnum , 'String' ) ) ;
mdsopen ( shotnum ) ;

try

g bolo=mdsdata ( '\results::btomo:emissivity' ) ;

comment=mdsdata ( '\results::btomo:comment' ) ;
delchans=mdsdata ( '\results::btomo:nchord' ) ;
Ptot=mdsdata ( '\results::btomo:prad_tot' ) ;
lambda=mdsdata ( '\results::btomo:lambda' ) ;
rpix=mdsdata ( '\results::btomo:rmesh' ) ;
zpix=mdsdata ( '\results::btomo:zmesh' ) ;
tvec=mdsdata ( '\results::btomo:time' ) ;
c h i s q=mdsdata ( '\results::btomo:chi_squared' ) ;
Pup=mdsdata ( '\results::btomo:prad_above' ) ;
Pbot=mdsdata ( '\results::btomo:prad_below' ) ;

catch

comment=' Data problem, sorry.'

end

mdsclose ;
i f length ( delchans )<20
delchans1=int2str ( delchans ’ ) ; ;
delchans2 = [ ] ;

else

delchans1=int2str ( delchans ( 1 : 1 9 ) ’ ) ;
delchans2 = int2str ( delchans ( 2 0 : length ( delchans ) ) ’ ) ;

APPENDIX B. THE PROGRAM CODE 81 end

nr=length ( rpix ) ; nz=length ( zpix ) ;
t s t e p s=length ( tvec ) ;
g bolo=reshape ( g bolo , nz , nr , t s t e p s ) ;
set (bm comment , 'String' , comment ) ;
pause (3)
set ( bm delchans , 'String' , delchans1 ) ;
set ( bm delchans2 , 'String' , delchans2 ) ;
set ( bm lambda , 'String' , num2str( lambda , 4 ) ) ;
i f t s t e p s==0

   set (bm comment , 'String' , 'No data, sorry.' ) ;

return ;

end

i f isnan (Pup ) Pup=NaN∗ ones ( size ( tvec ) ) ; end
i f isnan ( Pbot ) Pbot=NaN∗ ones ( size ( tvec ) ) ; end

set (bm comment , 'String' , '  ...wait please, making frames...' ) ;

maxlev=max( g bolo ( : ) ) ;
c l e v s =[0: wperc : maxlev ] ;
b f r =5;  % number of empty frames
for t i =1: t s t e p s

cla

[ con , hcon]=contour ( rpix , zpix , g bolo ( : , : , t i ) , c l e v s , 'k' ) ;
hold on ;
%  p l o t ( [ . 6 1 . 1 5 ] , [ lambda ( t i )/100 lambda ( t i ) / 1 0 0 ] , ’ b ’ ) ;
%  p l o t ( [ Ptot ( t i )/100 Ptot ( t i ) / 1 0 0 ] , [ −. 8 . 8 ] , ’ b ’ ) ;
%  p l o t ( [ . 6 1 . 1 5 ] , [ zmag( t i )/100 zmag( t i ) / 1 0 0 ] , ’m’ ) ;
%  p l o t ( [ rmag ( t i )/100 rmag ( t i ) / 1 0 0 ] , [ −. 8 . 8 ] , ’m’ ) ;
text ( . 6 5 , −. 6 8 , 'FRAME' ) ; text ( . 9 5 , −. 6 8 , 'TIME' ) ;
text ( . 7 8 , −. 7 5 , int2str ( t i ) ) ; text ( 1 , −. 7 5 , num2str( tvec ( t i ) , 4 ) ) ;
text ( . 6 5 , . 7 4 , [ int2str ( wperc /1000) ' kW.m^{-3}/cont.' ] ) ;
hold o f f ;
caxis ( [ 0 maxlev ] ) ;
axis equal ; shading f l a t ;
axis ( [ . 6 1 . 1 5 −. 8 . 8 ] ) ;
t i t l e ( [ 'shot #' int2str ( shotnum ) ] ) ;
xlabel ( 'r [cm]' ) ;
ylabel ( 'z [cm]' ) ;

set ( bm framenum , 'String' , num2str( t i ) ) ;
set ( bm Ptot , 'String' , num2str( Ptot ( t i ) / 1 0 0 0 , 4 ) ) ;
set (bm Pup , 'String' , num2str(Pup( t i ) / 1 0 0 0 , 4 ) ) ;
set ( bm Pbot , 'String' , num2str( Pbot ( t i ) / 1 0 0 0 , 4 ) ) ;
set ( bm chisq , 'String' , num2str( c h i s q ( t i ) , 4 ) ) ;

i f t i ==1 % moviein i s not necessary in MatLab6
moma=moviein ( t s t e p s+b f r , bm movieaxis ) ;

end

moma ( : , t i )=getframe ( bm movieaxis ) ;

end

cla ;
for blanc =1: b f r  % to d i s t i n g u i s h the end in movie r e p l y
moma ( : , t s t e p s+blanc)=getframe ( bm movieaxis ) ;

end

set ( bm framenum , 'String' , '1' ) ;

bolmov ( 'framenum' )

return ;

case 'play'
set (bm comment , 'String' , 'Playing movie...' ) ;

repeat=eval ( get ( bm repeat , 'String' ) ) ;
speed=eval ( get ( bm speed , 'String' ) ) ;
movie(moma, repeat , speed ) ;

bolmov ( 'framenum' )

return

case 'chan>'

framenum=eval ( get ( bm framenum , 'String' ))+1;
set ( bm framenum , 'String' , int2str ( framenum ) )

bolmov ( 'framenum' ) ;

return

case 'chan<'

framenum=eval ( get ( bm framenum , 'String' ))−1;
set ( bm framenum , 'String' , int2str ( framenum ) )

bolmov ( 'framenum' ) ;

return

case 'framenum'
set (bm comment , 'String' , '........' ) ;

t i=eval ( get ( bm framenum , 'String' ) ) ;

i f t i <1 t i =1;

set (bm comment , 'String' , '..hey, don''t exaggerate!' ) ; pause ( 1 ) ;

end

i f t i>t s t e p s t i=t s t e p s ;

set (bm comment , 'String' , '..hey, don''t exaggerate!' ) ; pause ( 1 ) ;

end drawnow cla

contour ( rpix , zpix , g bolo ( : , : , t i ) , c l e v s , 'k' ) ;
hold on ;
%  p l o t ( [ . 6 1 . 1 5 ] , [ lambda ( t i )/100 lambda ( t i ) / 1 0 0 ] , ’ b ’ ) ;
%  p l o t ( [ Ptot ( t i )/100 Ptot ( t i ) / 1 0 0 ] , [ −. 8 . 8 ] , ’ b ’ ) ;
%  p l o t ( [ . 6 1 . 1 5 ] , [ zmag( t i )/100 zmag( t i ) / 1 0 0 ] , ’m’ ) ;
%  p l o t ( [ rmag ( t i )/100 rmag ( t i ) / 1 0 0 ] , [ −. 8 . 8 ] , ’m’ ) ;
text ( . 6 5 , −. 6 8 , 'FRAME' ) ; text ( . 9 5 , −. 6 8 , 'TIME' ) ;
text ( . 7 8 , −. 7 5 , int2str ( t i ) ) ; text ( 1 , −. 7 5 , num2str( tvec ( t i ) , 4 ) ) ;
text ( . 6 5 , . 7 4 , [ int2str ( wperc /1000) ' kW.m^{-3}/cont.' ] ) ;
hold o f f ;
caxis ( [ 0 maxlev ] ) ;
axis equal ; shading f l a t ;
axis ( [ . 6 1 . 1 5 −. 8 . 8 ] ) ;
t i t l e ( [ 'shot #' int2str ( shotnum ) ] ) ;
xlabel ( 'r [cm]' ) ;
ylabel ( 'z [cm]' ) ;
set ( bm framenum , 'String' , num2str( t i ) ) ;
set ( bm Ptot , 'String' , num2str( Ptot ( t i ) / 1 0 0 0 , 4 ) ) ;
set ( bm lambda , 'String' , num2str( lambda , 4 ) ) ;
set (bm Pup , 'String' , num2str(Pup( t i ) / 1 0 0 0 , 4 ) ) ;
set ( bm Pbot , 'String' , num2str( Pbot ( t i ) / 1 0 0 0 , 4 ) ) ;
set ( bm chisq , 'String' , num2str( c h i s q ( t i ) , 4 ) ) ;
set (bm comment , 'String' , 'Ready.' ) ;

return end

B.3.2 fabreader.m

function fabreader ( shot ) ;

disp ( 'FABCAT MDS Reader. Release 1.4' ) ;

i f nargin < 1
s t a r t s h o t = 19000;
endshot = 21739;

APPENDIX B. THE PROGRAM CODE 84 else

s t a r t s h o t = shot ;
endshot = shot ;

end

index = 1;
for shot = s t a r t s h o t : endshot

try

shotnumber ( index ) = shot ;
exist ( index ) = 0 ;
xpoint ( index ) = NaN;
s t a t u s ( index ) = NaN;
lambda ( index ) = NaN;
chi2 ( index ) = NaN;
mdsopen ( 'results' , shot ) ;
disp ( [ 'Analyzing shot ' , num2str( shot ) ] ) ;

method = mdsdata ( '\results::btomo:method' ) ;

disp ( method ) ;
i f method (1:17) == ( 'FBTR with x-point' )
xpoint ( index ) = 1 ;

else

xpoint ( index ) = 0 ;

end

lambda ( index ) = mdsdata ( '\results::btomo:lambda' ) ;
chim = mdsdata ( '\results::btomo:chi_squared' ) ;

chi2 ( index ) = mean( chim ) ;
i f lambda ( index ) < 0 . 0 1 | abs ( chi2 ( index ) − 1) > 0.05
s t a t u s ( index ) = 0 ;

else

s t a t u s ( index ) = 1 ;

end

exist ( index ) = 1 ;

r e s u l t = mdsdata ( '\results::btomo:comment' ) ;

disp ( r e s u l t ) ;

catch end

index = index + 1;

mdsclose end

save f a b r e s u l t s shotnumber exist xpoint s t a t u s lambda chi2

B.3.3 fabanalyze.m

disp ( 'FABCAT Result Analyzer Release 1.4' ) ;

try

load r e s u l t s

catch

   disp ( 'File fabresults.mat not found. Run fabreader.m first.' ) ;

end

disp ( [ 'Shots out of interval ' , num2str( shotnumber ( 1 ) ) , …
' to ' , num2str( shotnumber ( length ( shotnumber ) ) ) ] ) ;
disp ( [ 'There are ' , num2str(sum( exist ) ) , ' shots available' ] ) ;
plus = 0;
for i = 1 : length ( xpoint )
i f ˜ isnan ( xpoint ( i ) )
i f xpoint ( i ) == 1
plus = plus + 1;

end end

end

disp ( [ 'There are ' , num2str( plus ) , ' shots with xpoint available' ] ) ;

conv = 0;
for i = 1 : length ( s t a t u s )
i f ˜ isnan ( s t a t u s ( i ) )
i f s t a t u s ( i ) == 1
conv = conv + 1;

end APPENDIX B. THE PROGRAM CODE 86

end end

disp ( [ 'There are ' , num2str(conv ) , ' shots with good reconstruction available' ] ) ;

B.3.4 fabplot.m

function f a b p l o t ( shot ) ;

disp ( 'FABCAT radiation Plotter. Release 1.4' ) ;

mdsopen ( 'results' , shot ) ;
disp ( [ 'Analyzing shot ' , num2str( shot ) ] ) ;

Pxpoint = mdsdata ( '\results::btomo:prad_above' ) ;
Pbelow = mdsdata ( '\results::btomo:prad_below' ) ;
Ptot = mdsdata ( '\results::btomo:prad_tot' ) ;
time = mdsdata ( '\results::btomo:time' ) ;

mdsclose ;
mdsopen ( shot ) ;
Ip=mdsdata ( '\results::I_p' ) ;  % plasma current

 tIp = mdsdata ( 'dim_of(\results::I_p)' ) ;
 uIp = mdsdata ( 'units_of(\results::I_p)' ) ;

nel = mdsdata ( '\results::fir:n_average' ) ;  % plasma average density

 t n e l = mdsdata ( 'dim_of(\results::fir:n_average)' ) ;
 unel = mdsdata ( 'units_of(\results::fir:n_average)' ) ;

[ tPoh , Poh]= tcvget ( 'POHM' ) ;

figure ( 'Name' , 'Radiation of x-point' ) ;

plot ( Pxpoint ) ;

figure ( 'Name' , 'Radiation below the x-point' ) ;

plot ( Pbelow ) ;

figure ( 'Name' , 'Total radiation' ) ;

plot ( Ptot ) ;

figure ( 'Name' , 'Plasma current' ) ;

plot ( Ip ) ;

figure ( 'Name' , 'Plasma density' ) ;

plot ( nel ) ;

figure ( 'Name' , 'Ohmic power' ) ;

plot (Poh ) ;

figure ( 'Name' , 'Plasma core flux surface' ) ;

tcvview ( 'pvt' , shot , . 5 )

figure ( 'Name' , 'Comparison between Ptot and Pohm' ) ;

plot ( tPoh , Poh , 'g' , time , Ptot , 'b' ) ;

legend ( 'Ohmic power' , 'Total radiated power' ) ;

axis ( [ 0 time ( length ( time ) ) + 0 . 1 0 max(max( Ptot ) ,max(Poh ) ) ] ) ;

figure ( 'Name' , 'Comparison between radiated powers' ) ;

plot ( time , Ptot , 'r' , time ,8∗ Pbelow , 'b' , time ,8∗ Pxpoint , 'g' ) ;

legend ( 'Total radiated power' , 'Power radiated below x-point' , 'xpoint radiation' ) ;

axis ( [ 0 time ( length ( time ) ) + 0 . 1 0 max( Ptot ) ] ) ;
mdsclose ;