BLCoil.hh

Go to the documentation of this file.

//      BLCoil.hh
/*
This source file is part of G4beamline, http://g4beamline.muonsinc.com
Copyright (C) 2003,2004,2005,2006 by Tom Roberts, all rights reserved.

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

http://www.gnu.org/copyleft/gpl.html
*/

#ifndef BLCOIL_HH
#define BLCOIL_HH

#include <map>

#include "globals.hh"
#include "BLCommand.hh"
#ifdef STUB
#include "NTuple.hh"
#endif // STUB

/**     BLCoil creates and implements a coil (part of a solenoid)
 *
 *      A BLCoil consists only of the parameters describing the coil geometry
 *      and a field map for unit current density. It is not a BLElement, and
 *      cannot be placed geometrically -- see BLSolenoid for that.
 *
 *      Cylindrical symmetry around z is assumed. The field is approximated
 *      by computing the field from nSheets infinitesimally-thin current
 *      sheets interpolating from innerRadius to OuterRadius, each length
 *      long. The solution to Maxwell's Equations for each sheet is used.
 *
 *      The field has been verified to be correct within 0.2% for both a
 *      long solenoid and for a pair of Helmholtz coils.
 *
 *      Effect of nSheets on accuracy: consider the coil innerRadius=330.0*mm,
 *      outerRadius=505.0*mm, length=167.0*mm. Compute the maximum error
 *      for 1000 random points, excluding the region within 2 cm of the coils;
 *      use 99 sheets as the "correct" value. In the region 0<=r<2*outerRadius,
 *      0<z<10*outerRadius, for 0.1% accuracy one needs 11 sheets, and for 1%
 *      one needs 4 sheets. In the region 0<=r<innerRadius/2, same z, for 0.1%
 *      accuracy one needs 5 sheets, and for 1% one needs 2 sheets. Note that
 *      inside the coils this approximation of using thin sheets is not very
 *      good; but interpolating on a grid and coordinating sheet positions
 *      with the grid is probably about as good as one can do -- at least the
 *      values will be smooth (which is not true for the bare computation --
 *      there's a singularity at every sheet, and Bz jumps from +infinity
 *      to -infinity across the sheet). 
 *
 *      Note that Bz does not decrease to 0.1% of its central value until
 *      z>10*(innerRadius+outerRadius)/2; it decreases to 1% at ~5 times that
 *      average radius.
 *
 *
 *      BLCoil can also read a mapFile defining the field (rather than doing
 *      the above computation). For instance, it could be from opera-2d.
 *      The mapFile gives tables of Bz and Br on a grid in the z-r plane.
 *      mapFile format:
 *              lines beginning with # are comments, and are ignored.
 *              The first line contains the parameters in name=value form:
 *                      z0      initial value of z (mm)
 *                      dz      grid interval along z (mm)
 *                      nz      number of grid points along z
 *                      dr      grid interval along r (mm)
 *                      nr      number of grid points along r
 *                      reflectZ set nonzero to reflect in the z=0 plane
 *                              (inverting Br)
 *                      current the current for this field (A/mm^2)
 *              All parameters are required. If reflectZ is nonzero then
 *              z0 must be 0.0 (only z>=0 is contained in the file).
 *              The point with local coordinates (z=0,r=0) is what gets placed.
 *
 *              Then follow the two maps in any order. Each map begins
 *              with a line giving its name (Bz or Br) in column 1, followed
 *              by nz lines each containing nr floats (Tesla). The first
 *              value in each line is for r=0. Values are separated by one
 *              or more spaces (the first value on each line may but need not
 *              be preceeded by spaces).
 *      The main reason to use a mapFile is to accurately include the effects
 *      of iron. That implies that the map may not be accurate for values
 *      of current other than that for which the map was generated.
 *
 *      See BLCMDfieldmap for a more general EM field read from a map-file.
 *
 *
 *      BLCoil was separated from BLSolenoid so multiple solenoids with the
 *      same geometry (coil) could share a single (large) field map for the
 *      coil -- the B field is scaled by current/currentOfMap, where 
 *      currentOfMap is 1.0 for the computed field, and is the current 
 *      parameter from a mapFile.
 *
 *      This class does not use BLFieldMap becuase it predates it, and
 *      because the map must include the parameters so construct() will
 *      know whether or not the map corresponds to the current values.
 *      You can, of course, build a solenoid out of a BLTubs and a
 *      BLFieldMap.
 **/
class BLCoil {
        static std::map<G4String,BLCoil*> mapCoil;
        G4String name;
        G4double innerRadius;
        G4double outerRadius;
        G4double length;
        G4int nSheets;
        G4String material;
        G4double tolerance;
        G4double maxR;
        G4double minZ;
        G4double maxZ;
        G4double dR;
        G4double dZ;
        G4int nR;
        G4int nZ;
        G4String filename;
        G4String mapFile;
        G4int exactComputation;
        float *mapBr;
        float *mapBz;
        G4double sheetR0;
        G4double sheetDR;
        G4double norm;
        G4bool goodMap;
        G4bool reflectZ;
        G4double currentOfMap;
        friend class BLCMDcoil;
        friend class BLCMDsolenoid;
        friend class SolenoidField;
public:
        /// Default constructor.
        BLCoil();

        /// Complex Constructor. All arguments that default to 0 are related
        /// to generating the field map; those that are 0 will be determined
        /// by tolerance. To construct the map the minimum info is tolerance
        /// and maxR; if either of them is 0, no map will be generated or used.
        /// NOTE: without a map, field values inside the coils will fluctuate
        /// wildly, because points on the sheets are singular; map construction
        /// is designed to minimize this. Without a map, nSheets must be given.
        /// tolerance is a fraction of the best value at r=0,z=0; values
        /// smaller than about 0.001 are unrealistic (and expensive!).
        BLCoil(G4String _name, G4double _innerRadius, G4double _outerRadius,
                G4double _length, G4int _nSheets=0, G4double _tolerance=0.0,
                G4double _maxR=0.0, G4double _maxZ=0.0, G4double _dR=0.0,
                G4double _dZ=0.0, G4int _nR=0, G4int _nZ=0, 
                G4String _filename="") {
                        name=_name; innerRadius=_innerRadius;
                        outerRadius=_outerRadius; length=_length;
                        nSheets=_nSheets; tolerance=_tolerance; maxR=_maxR;
                        maxZ=_maxZ; dR=_dR; dZ=_dZ; nR=_nR; nZ=_nZ; 
                        filename = _filename;
                        exactComputation = 0;
                        mapBr=0; mapBz=0; sheetR0=0.0; sheetDR=0.0; norm=0.0; 
                        goodMap=false; 
                        if(nSheets > 0) { 
                                sheetDR=(outerRadius-innerRadius)/nSheets;
                                sheetR0=innerRadius + sheetDR/2.0;
                        }
                        if(filename == "") filename = name + ".dat";
                }

        /// Destructor.
        virtual ~BLCoil();

        /// Copy constructor.
        BLCoil(const BLCoil& r);

        /// printCoil() will print a description
        void printCoil();

        /// find() returns a pointer to the named BLCoil.
        static BLCoil *find(G4String name);

        /// addField() adds the field value for the given point[] into field[].
        /// Note that point[] is relative to the center of the coil.
        /// Units for current are Amp/mm^2.
        void addField(const G4double point[4], G4double field[6],
                                        G4double currentAmpPerMm2) const;

        /// generateFieldMap() will generate the field map (for execution
        /// speedup). If just maxR and tolerance were given, all other
        /// map parameters will be automatically determined so that the
        /// estimated error is less than the specified tolerance.
        /// The field map is cached in a disk file which is automatically
        /// recreated if anything changes (may take several minutes).
        void generateFieldMap();

        /// readMapFile() will read the mapFile.
        void readMapFile();

        /// displayMapValues() will display the 4 map values around a point.
        void displayMapValues(G4double r, G4double z);
private:
        /// getUnitField() will return Br and Bz for unit current.
        void getUnitField(G4double r, G4double z, G4double& Br, G4double& Bz);

        /// determineMaxR() will determine maxR based on tolerance.
        /// Uses nSheets=99 to determine center field, and walks out
        /// far enough so (Bz/Bz0)<0.5*tolerance.
        void determineMaxR();

        /// determineMaxZ() will determine maxZ based on tolerance.
        /// Uses nSheets=99 to determine on-axis field, and walks out
        /// far enough so (Bz/Bz0)<0.5*tolerance.
        void determineMaxZ();

        /// determineNsheets() will determine nSheets, given maxZ and 
        /// tolerance. Increases nSheets until maxError < tolerance/4
        /// (the other 3/4 of tolerance is for mapping errors). Checks
        /// just a handful of points selected to be representative.
        void determineNsheets();

        /// estimateMapError() will estimate the maximum error in mapping,
        /// given maxZ, nSheets, dR, dZ, nR, nZ, and an allocated but
        /// uninitialized fieldMap. It generates 1000 random points and
        /// returns the maximum error; quits early if tolerance is
        /// exceeded. Error is (currentMap - best)/best(r=0,z=0), where 
        /// best is computed with no map and nSheets=99. For each point
        /// it generates the 4 surrounding values of fieldMap (this is faster
        /// than generating the entire map first).
        /// Returns the maximum error of either Bz or Br, normalized to
        /// Bz(r=0,z=0,nSheets=99).
        G4double estimateMapError();

        /// computeError() computes the error for a single point.
        /// Error is (currentMap - best)/best(r=0,z=0), where best should be
        /// an identical BLCoil with no map and nSheets=99. If makeMap is true,
        /// it generates the 4 surrounding values of the fieldMap.
        G4double computeError(G4double r, G4double z, BLCoil &best,
                                                        G4bool makeMap);

        /// getSheetField() returns the B field from an infinitesimally-thin
        /// sheet of current. Arguments:
        /// r and z are the cylindrical coordinates of the point in question
        /// relative to center of the coil (axis along z), mm.
        /// Br and Bz return the computed B field at that point (geant4 units)
        /// sheetR is the radius of the current sheet (mm)
        /// sheetLen is the length of the current sheet (mm)
        /// currentAmpPerMm is the current of the sheet (Amp per linear mm)
        void getSheetField(G4double r, G4double z, G4double& Br, G4double& Bz,
                G4double sheetR, G4double sheetLen, G4double currentAmpPerMm)
                                                                        const;
};

#endif // BLCOIL_HH