// AK20100820 Modified from G4beamline v2.03 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 MYBLCOIL_HH
#define MYBLCOIL_HH

#include <map>

#include "IBLFuncs.hxx"

namespace COMET {

/**	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 IBLCoil {
    static std::map<std::string,IBLCoil*> mapCoil;
    std::string name;
    Double_t innerRadius;
    Double_t outerRadius;
    Double_t length;
    Int_t nSheets;
    std::string material;
    Double_t tolerance;
    Double_t maxR;
    Double_t minZ;
    Double_t maxZ;
    Double_t dR;
    Double_t dZ;
    Int_t nR;
    Int_t nZ;
    std::string filename;
    std::string mapFile;
    float *mapBr;
    float *mapBz;
    Double_t sheetR0;
    Double_t sheetDR;
    Double_t norm;
    Bool_t goodMap;
    Bool_t reflectZ;
    Double_t currentOfMap;
    public:
        /// Default constructor.
        IBLCoil();

        /// 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!).
        IBLCoil(std::string _name, Double_t _innerRadius, Double_t _outerRadius,
            Double_t _length, Int_t _nSheets=0, Double_t _tolerance=0.0,
            Double_t _maxR=0.0, Double_t _maxZ=0.0, Double_t _dR=0.0,
            Double_t _dZ=0.0, Int_t _nR=0, Int_t _nZ=0,
        std::string _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;
            mapBr=0;
            mapBz=0;
            sheetR0=0.0;
            sheetDR=0.0;
            norm=0.0;
            currentOfMap=1.0;

            goodMap=false;
            if(nSheets > 0) {
                sheetDR=(outerRadius-innerRadius)/nSheets;
                sheetR0=innerRadius + sheetDR/2.0;
            }
            if(filename == "") filename = name + ".dat";
        }

        /// Destructor.
        virtual ~IBLCoil();

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

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

        /// find() returns a pointer to the named IBLCoil.
        static IBLCoil *find(std::string name);

        /// find() returns a pointer to the named IBLCoil.
        static void registerCoil(IBLCoil* coil){
            mapCoil[coil->name]=coil;
        }

        /// 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 Double_t point[4], Double_t field[6],
            Double_t 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(Double_t r, Double_t z);

        /// Get the length of this coil
        double GetLength()const{return length;}
        /// Get the inner radius of this coil
        double GetInnerRadius()const{return innerRadius;}
        /// Get the outer radius of this coil
        double GetOuterRadius()const{return outerRadius;}
        /// Get the maximum radius to use this coil for field calculations
        double GetMaxR()const{return maxR;}
        /// Get the maximum Z to use this coil for field calculations
        double GetMaxZ()const{return maxZ;}

    private:
        /// getUnitField() will return Br and Bz for unit current.
        void getUnitField(Double_t r, Double_t z, Double_t& Br, Double_t& 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).
        Double_t estimateMapError();

        /// computeError() computes the error for a single point.
        /// Error is (currentMap - best)/best(r=0,z=0), where best should be
        /// an identical IBLCoil with no map and nSheets=99. If makeMap is true,
        /// it generates the 4 surrounding values of the fieldMap.
        Double_t computeError(Double_t r, Double_t z, IBLCoil &best,
            Bool_t 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(Double_t r, Double_t z, Double_t& Br, Double_t& Bz,
            Double_t sheetR, Double_t sheetLen, Double_t currentAmpPerMm)
            const;
};
}
#endif                                                // MYBLCOIL_HH