superarray.h

Go to the documentation of this file.

#pragma once
#include <stdlib.h>
#include <fstream>
#include "vector.h"
#include "resutil.h"
#include "objectshape.h"
#include "error.h"
#include <vector>
#include "gridentry.h"
 
// #define INDEX_TYPE long long int
namespace GOAT
{
   namespace raytracing 
   {
       typedef long long int INDEX_TYPE; 
      template <class T> class SuperArray
      {
         public :
         SuperArray ();
         SuperArray(double r0, INDEX_TYPE nx, INDEX_TYPE ny, INDEX_TYPE nz, const int typ=IN_OBJECT);
         SuperArray(double r0, INDEX_TYPE nx, INDEX_TYPE ny, INDEX_TYPE nz, std::vector<ObjectShape*>Obj, int numObjs, const bool isAbsolute=false, const int typ=IN_OBJECT);
         SuperArray(const SuperArray& S)
         {
             Error = NO_ERRORS;
             
             type = S.type;
             ywerte = S.ywerte;
             zwerte = S.zwerte;
             numObjs = S.numObjs;
             Obj = S.Obj;
             G = S.G;
             K = S.K;
             Pul = S.Pul;
             n = S.n;
             nges = S.nges;
             d = S.d;
             r0 = S.r0;
             isequal = S.isequal;
             iscleared = S.iscleared;
             pc = S.pc;
             H = S.H;
             R = S.R;
         }
         ~SuperArray () { clear();}
         
         void copy (const SuperArray &S); 
         void add (const SuperArray& S);
         void sub (const SuperArray& S);
      //   friend double sumabs (const SuperArray &S); ///< The absolute values of the cell contents are summed up
      //   friend double sumabs2 (const SuperArray &S); ///< The squared absolute values of the cell contents are summed up
      //   friend double abs2sum (const SuperArray &S);
         bool addInc(std::vector<ObjectShape*>Obj,int numObj, const bool isAbsolute=false); // numObjs Einschluesse hinzufuegen
                                           // isAbsolute=true => Orts-/Groessenangaben in absoluten Koordinaten 
         bool addInc (ObjectShape *E,const bool isAbsolute=false);  // Einschluss hinzufuegen (isAbsolute s.o.) 
         void removeObject(int i); 
         void removeAllObjects();
         void setActive(int i);
         void setInactive(int i);
         maths::Vector<INDEX_TYPE> gitterpunkt (maths::Vector<double> P);
         bool inObject (maths::Vector<double> P, int i); 
         bool inObject (maths::Vector<INDEX_TYPE> Pi, int i); 
         bool inObject (INDEX_TYPE ix, INDEX_TYPE iy, INDEX_TYPE iz, int i); 
         T& operator () (INDEX_TYPE ix, INDEX_TYPE iy, INDEX_TYPE iz); 
         T& operator () (int i, INDEX_TYPE ix, INDEX_TYPE iy, INDEX_TYPE iz, bool isObjectCoordinate=true); 
         T& operator () (maths::Vector<INDEX_TYPE> Pi); 
         T& operator () (int i,maths::Vector<INDEX_TYPE> Pi); 
         T& operator () (maths::Vector<double> P); 
         T& operator () (int i, maths::Vector<double> P); 
         
          void makeReal ();
         void fill(const T &x); 
         SuperArray& operator = (const SuperArray &S); 
         void clear(); 
         void reinit(double r0, INDEX_TYPE n);
         void reinit(double r0, INDEX_TYPE nx, INDEX_TYPE ny, INDEX_TYPE nz);
     
         void allockugel(); 

         maths::Vector<INDEX_TYPE> kugelindex(maths::Vector<INDEX_TYPE> Pi); 
         T kugelwert(maths::Vector<int> Pi); 
         T kugelwert(int ix, int iy, int iz); 
 
         int Error;  
         std::vector<ObjectShape*> Obj; 
         int numObjs; 
         int type; 
         std::vector<int> ywerte;
         std::vector<std::vector<int> > zwerte;
         std::vector <std::vector <std::vector <std::vector <T> > > > G; 
         std::vector < std::vector < std::vector <T> > >  K;
         T dummy;
         std::vector<maths::Vector<INDEX_TYPE>> Pul;
         std::vector<maths::Vector<INDEX_TYPE>> n; 
         maths::Vector<INDEX_TYPE> nges; 
 
         maths::Vector<double> d;  
         double r0; 
         bool isequal;
         bool iscleared = true; 
         T pc;        
         maths::Matrix<double> H, R; 
         bool write(std::string fname);
         bool read(std::string fname);

         void setNumberOfCellsPerDirection(INDEX_TYPE no) {
             // Ensure the parameter is valid and reinitialize the grid with the given number of cells per direction  
             if (no > 0) {
                 reinit(r0, no, no, no);
             }
             else {
                 Error = INVALID_PARAMETER; // Set an error code for invalid input  
             }
         }
 
         
             void arrangeGrids();
    };
 
    /* Specialization of Superarray for use together with std::vector<gridEntry> */
    template<> void SuperArray<std::vector<GOAT::raytracing::gridEntry> >::clear();
    template <class T> T dummy;
    
 
    /*------------------------- IMPLEMENTATION --------------------------------------*/
 
    template<> bool SuperArray<GOAT::maths::Vector<std::complex<double> > >::write(std::string fname);
 
    template <class T> SuperArray<T>::SuperArray()
    {
        H = maths::unity();
        R = maths::unity();
        Error = NO_ERRORS;
        numObjs = 0;
        isequal = false;
 
        type = IN_HOST;        
    }
 
    template <class T> SuperArray<T>::SuperArray(double r0, INDEX_TYPE nx, INDEX_TYPE ny, INDEX_TYPE nz, const int typ)
    {
        Error = NO_ERRORS;
        double b = 2.0 * r0;
        isequal = false;
        numObjs = 0;
        this->r0 = r0;
 
 
 
        nges = maths::Vector<INDEX_TYPE>(nx, ny, nz);
        d = maths::Vector<double>(b / (double)nx, b / (double)ny, b / (double)nz);
        iscleared = true;
        type = typ;
        if (type & IN_HOST)
            allockugel();
    }
 
    template <class T> SuperArray<T>::SuperArray(double r0, INDEX_TYPE nx, INDEX_TYPE ny, INDEX_TYPE nz, std::vector<ObjectShape*> Obj, int numObjs, const bool isAbsolute, const int typ)
    {
        Error = NO_ERRORS;
        double b = 2.0 * r0;
        isequal = isAbsolute;
        this->numObjs = 0;
        this->r0 = r0;
        nges = maths::Vector<INDEX_TYPE>(nx, ny, nz);
        type = typ;
        d = maths::Vector<double>(b / nx, b / ny, b / nz);
        iscleared = true;
        if (type & IN_HOST)
        {
            allockugel();
        }
        addInc(Obj, numObjs, isAbsolute);
    }
 
 
    template <class T> bool SuperArray<T>::addInc(std::vector<ObjectShape*> Obj, int numObjs, const bool isAbsolute)
    {
        bool ok = true;
        for (int i = 0; (i < numObjs)/* && (canCalc)*/ && (ok); i++) ok = addInc(Obj[i], isAbsolute);
        return ok;
    }
 
    template<class T>
    inline void SuperArray<T>::setInactive(int i)
    {
        removeObject(i);
    }
 
    template<class T>
    inline void GOAT::raytracing::SuperArray<T>::setActive(int i)
    {
        if (G[i].size() == 0) G[i].resize(n[i][0] + 1);
            for (INDEX_TYPE ix = 0; ix < n[i][0] + 1; ix++)
            {
                G[i][ix].resize(n[i][1] + 1);
 
                for (INDEX_TYPE iy = 0; iy < n[i][1] + 1; iy++)
                {
                    G[i][ix][iy].resize(n[i][2] + 1);
                }
            }
    }
 
    template <class T> maths::Vector<INDEX_TYPE> SuperArray<T>::gitterpunkt(maths::Vector<double> P)
    {
        maths::Vector<INDEX_TYPE> pi, ph;
        maths::Vector<double> h;
        h = H * P + maths::Vector<double>(r0, r0, r0);
        pi = maths::Vector<INDEX_TYPE>((INDEX_TYPE)floor(h[0] / d[0]), (INDEX_TYPE)floor(h[1] / d[1]), (INDEX_TYPE)floor(h[2] / d[2]));
        /*
        if (type & IN_HOST) //????????
        {
            pi = ph;
        }
        else
        {
            pi = ph;
        }*/
        return pi;
    }

    template <class T> void SuperArray<T>::arrangeGrids()
    {
        clear();
        maths::Vector<INDEX_TYPE> pul, por, N, hn;
        maths::Vector<double> h, O;
        O = maths::Vector<double>(-r0, -r0, -r0);
        d = 2.0 * r0 * maths::Vector<double>(1.0 / (double)nges[0], 1.0 / (double)nges[1], 1.0 / (double)nges[2]);
        std::cerr << "[arrangeGrids] d=" << d << "\tr0=" << r0 << "\tnges=(" << nges[0] << "," << nges[1] << "," << nges[2] << ")" << std::endl;
        
        for (int i = 0; i < numObjs; i++)
        {
            Obj[i]->initQuad();
;            h = ceil(ediv(Obj[i]->por, d)) - floor(ediv(Obj[i]->pul, d));
            hn = maths::Vector<INDEX_TYPE>((INDEX_TYPE)h[0] + 1, (INDEX_TYPE)h[1] + 1, (INDEX_TYPE)h[2] + 1); // Gr��e des 3D-Gitters in die drei Koordinatenrichtungen
            n.push_back(hn);
            
            h = floor(ediv(Obj[i]->pul + maths::Vector<double>(r0, r0, r0), d));
            Pul.push_back(maths::Vector<INDEX_TYPE>((INDEX_TYPE)h[0], (INDEX_TYPE)h[1], (INDEX_TYPE)h[2]));
 
            if (Obj[i]->isActive())  // Ist der Einschluss �berhaupt inelastisch aktiv ? 
            {
                G[i].resize(n[i][0] + 1);
                for (INDEX_TYPE ix = 0; ix < n[i][0] + 1; ix++)
                {
                    G[i][ix].resize(n[i][1] + 1);
 
                    for (INDEX_TYPE iy = 0; iy < n[i][1] + 1; iy++)
                    {
                        G[i][ix][iy].resize(n[i][2] + 1);
                    }
                }
            }
        }
    }
    
    template <class T> bool SuperArray<T>::addInc(ObjectShape* E, const bool isAbsolute)
    {
       
        //SysMemInfo smi;
        //MemInfo mi;
        long int allocMem;
        maths::Vector<INDEX_TYPE> pul, por, N, hn;
        maths::Vector<double> h, O;
        O = maths::Vector<double>(-r0, -r0, -r0);
        double b = 2.0 * r0;
        d = maths::Vector<double>(b / (double)(nges[0] - 1), b / (double)(nges[1] - 1), b / (double)(nges[2] - 1));
      //   h = ceil(ediv(E->por, d)) - floor(ediv(E->pul, d));
        h = ceil(ediv(E->por-E->pul, d)) ;
 
        hn = maths::Vector<INDEX_TYPE>((INDEX_TYPE)h[0]+1, (INDEX_TYPE)h[1]+1, (INDEX_TYPE)h[2]+1); // Gr��e des 3D-Gitters in die drei Koordinatenrichtungen
 
        /* Berechne den tats�chlichen Bedarf */
        allocMem = sizeof(T***) + 2 * sizeof(maths::Vector<int>) + sizeof(ObjectShape*)
            + hn[0] * sizeof(T**)
            + hn[0] * hn[1] * sizeof(T*)
            + hn[0] * hn[1] * hn[2] * sizeof(T);
 
        {
            if (numObjs < 1)  // Es ist der erste Einschluss der hinzugef�gt wird
            {   
                G.resize(1);                        
                G.shrink_to_fit();
            }
            else
            {
                G.resize(numObjs + 1);
                G.shrink_to_fit();
            }
        }
      
        n.push_back(hn);
        Obj.push_back(E);
        h = floor(ediv(Obj[numObjs]->pul + maths::Vector<double>(r0, r0, r0), d));
        Pul.push_back(maths::Vector<INDEX_TYPE>((INDEX_TYPE)h[0], (INDEX_TYPE)h[1], (INDEX_TYPE)h[2]));
        
        if (E->isActive())  // Ist der Einschluss �berhaupt inelastisch aktiv ? 
        {
//            std::cout << "n[" << numObjs << "]=" << n[numObjs] << std::endl;
            G[numObjs].resize(n[numObjs][0] + 1);
            
            for (INDEX_TYPE ix = 0; ix < n[numObjs][0] + 1; ix++)
            {
                G[numObjs][ix].resize(n[numObjs][1] + 1);                
                
                for (INDEX_TYPE iy = 0; iy < n[numObjs][1] + 1; iy++)
                {
                    G[numObjs][ix][iy].resize(n[numObjs][2] + 1) ;                                       
                }
            }
        }
        
 
        int i = 0;
        numObjs++;
        iscleared = false;
        return true;
    }
 
    template<class T>
    inline void SuperArray<T>::removeObject(int i)
    {
 
        if (G[i].size() > 0)
        {
            for (int ix = n[i][0]; ix >= 0; ix--)
            {
                for (int iy = n[i][1]; iy >= 0; iy--)
                    G[i][ix][iy].clear();
                G[i][ix].clear();
            } // for ix    
            G[i].clear();
        } // if (G[i]!=NULL)
    }
    template<class T>
    inline void SuperArray<T>::removeAllObjects()
    {
        for (int i = 0; i < Obj.size(); i++)
        {
            removeObject(i);
        }
 
        Obj.clear();
        Obj.shrink_to_fit();
        numObjs = 0;
 
       
    }
 
 
    template <class T> bool SuperArray<T>::inObject(maths::Vector<double> P, int i) // da muss noch was gemacht werden
    {
        return Obj[i]->isInside(P);
    }
 
    template <class T> bool SuperArray<T>::inObject(maths::Vector<INDEX_TYPE> Pi, int i)
    {
        return (Obj[i]->isInside(emult(Pi, d)));
    }
 
    template <class T> bool SuperArray<T>::inObject(INDEX_TYPE ix, INDEX_TYPE iy, INDEX_TYPE iz, int i)
    {
        return (Obj[i]->isInside(emult(maths::Vector<double>(ix, iy, iz), d)));
    }
 
    template <class T> T& SuperArray<T>::operator () (INDEX_TYPE ix, INDEX_TYPE iy, INDEX_TYPE iz)
    {
        maths::Vector<INDEX_TYPE> Pi = maths::Vector<INDEX_TYPE>(ix, iy, iz);
        INDEX_TYPE i = 0;
        bool found = false;
        if (type == IN_OBJECT)
        {
            do
            {
                found = inObject(ix, iy, iz, i);
                i++;
            } while ((i < numObjs) && (!found));
            i--;
            if (!found)
            {
                Error = NOT_FOUND;
                return dummy;
            }
            else
            {
                Pi = Pi - Pul[i];
                Pi = kugelindex(Pi);
                if (Error == NO_ERRORS)
                    return G[i][Pi[0]][Pi[1]][Pi[2]];
                else
                {
                    Error = NOT_FOUND;
                    return dummy;
                }
            }
        }
        else
        {
            Pi = kugelindex(Pi);
            if (Error == NO_ERRORS)
                return K[Pi[0]][Pi[1]][Pi[2]];
            else
            {
                Error = NOT_FOUND;
                return dummy;
            }
        }
    }
 
    template <class T> T& SuperArray<T>::operator () (maths::Vector<INDEX_TYPE> Pi)
    {
        // dummy = maths::Vector<std::complex<double> >(0.0, 0.0, 0.0);
       
        if (type == IN_OBJECT)
        {
            int i = 0;
            bool found = false;
            do
            {
                found = inObject(Pi, i);
                i++;
            } while ((i < numObjs) && (!found));
            i--;
            
            if (!found)
            {
                Error = NOT_FOUND;
                return dummy;
            }
            else
            {
                Pi = Pi - Pul[i];
                Pi = kugelindex(Pi);                
                if (Error == NO_ERRORS)
                    return G[i][Pi[0]][Pi[1]][Pi[2]];
                else
                {
                    Error = NOT_FOUND;
                    return dummy;
                }
            }
        }
        else
        {
            Pi = kugelindex(Pi);
            if (Error == NO_ERRORS)
                return K[Pi[0]][Pi[1]][Pi[2]];
            else
            {
                Error = NOT_FOUND;
                return dummy;
            }
        }
    }
 
 
    template <class T> T& SuperArray<T>::operator () (int i, INDEX_TYPE ix, INDEX_TYPE iy, INDEX_TYPE iz, bool isEinKoord)
    {
        maths::Vector<INDEX_TYPE> Pi = maths::Vector<INDEX_TYPE>(ix, iy, iz);
        if (type == IN_OBJECT)
        {
            if (G[i].size() == 0)
            {
                Error = NOT_FOUND;
                return dummy;
            }
            if (!isEinKoord)  Pi = Pi - Pul[i];
            return G[i][Pi[0]][Pi[1]][Pi[2]];
        }
        else
        {
            Pi = kugelindex(Pi);
            if (Error == NO_ERRORS)
                return K[Pi[0]][Pi[1]][Pi[2]];
            else
            {
                Error = NOT_FOUND;
                return dummy;
            }
        }
    }
 
 
    template <class T> T& SuperArray<T>::operator () (int i, maths::Vector<INDEX_TYPE> Pi)
    {
        // T dummy;
        if (type == IN_OBJECT)
        {
            Pi = Pi - Pul[i];
            if (G[i].size() == 0) { Error = NOT_FOUND; return dummy; }
            if (Pi[0] < 0) { Error = NOT_FOUND; return dummy; }//maths::Vector<std::complex<double> > (0,0,0);
            if (Pi[1] < 0) { Error = NOT_FOUND; return dummy; } //maths::Vector<std::complex<double> > (0,0,0);
            if (Pi[2] < 0) {
     //           std::cout << "Pi[2]=" << Pi[2] << std::endl;
                  Error = NOT_FOUND; return dummy; } // maths::Vector<std::complex<double> > (0,0,0);
            if (Pi[0] >= n[i][0])
            {
       //         std::cout << "PROBLEM: Pi[0]=" << Pi[0] << "\tn[i][0]=" << n[i][0] <<  std::endl;
                Error = SUPERGITTER;
                 return dummy; //maths::Vector<std::complex<double> > (0,0,0);
            }
 
            if (Pi[1] >= n[i][1])
            {
      //          std::cout << "PROBLEM: Pi[1]=" << Pi[1] << "\tn[i][1]=" << n[i][1] <<  std::endl;
                Error = SUPERGITTER;
                 return dummy; //maths::Vector<std::complex<double> > (0,0,0);
            }
 
            if (Pi[2] >= n[i][2])
            {
        //       std::cout << "PROBLEM: Pi[1]=" << Pi[2] << "\tn[i][2]=" << n[i][2] <<  std::endl;
                Error = SUPERGITTER;
                 return dummy;
            }
 
            Error = NO_ERRORS;
           return G[i][Pi[0]][Pi[1]][Pi[2]];
        }
        else
        {
            Pi = kugelindex(Pi);
            if (Error == NO_ERRORS)
                return K[Pi[0]][Pi[1]][Pi[2]];
            else
            {
                Error = NOT_FOUND; return dummy;
            }
        }
    }
 
    template <class T> T& SuperArray<T>::operator () (maths::Vector<double> P)
    {
        int i;
        maths::Vector<int> Pi;
        maths::Vector<double> h;
        h = H * P - maths::Vector<double>(-r0, -r0, -r0);
        dummy = maths::Vector<std::complex<double> >(0.0, 0.0, 0.0);
        for (i = 0; i < 3; i++)
            Pi[i] = h[i] / d[i];
 
        if (type == IN_OBJECT)
        {
            bool found = false;
            i = 0;
            do
            {
                found = inObject(Pi, i);
                i++;
            } while ((i < numObjs) && (!found));
            i--;
            // if (!found) return maths::Vector<std::complex<double> > (INF,INF,INF);
            if (!found)
            {
                Error = NOT_FOUND;
                return dummy;
            }
            else
            {
                if (G[i] == NULL)
                {
                    Error = NOT_FOUND;
                    return dummy;
                }
                Pi = Pi - Pul[i];
                std::cout << "PI:" << Pi << std::endl;
                return G[i][Pi[0]][Pi[1]][Pi[2]];
            }
        }
        else
        {
            Pi = kugelindex(Pi);
            if (Error == NO_ERRORS)
                return K[Pi[0]][Pi[1]][Pi[2]];
            else
            {
                Error = NOT_FOUND;
                return dummy;
            }
        }
    }
 
    template <class T> T& SuperArray<T>::operator () (int i, maths::Vector<double> P)
    {
        maths::Vector<int> Pi;
        maths::Vector<double> h;
        h = H * P - maths::Vector<double>(-r0, -r0, -r0);
        T dummy;
        // dummy = maths::Vector<std::complex<double> >(0.0, 0.0, 0.0);
        Error = NO_ERRORS;
        for (int j = 0; j < 3; j++)
            Pi[j] = h[j] / d[j];
 
        if (type == IN_OBJECT)
        {
            if (G[i] == NULL)
            {
                Error = NOT_FOUND;
                return dummy;
            }
            Pi = Pi - Pul[i];
            return G[i][Pi[0]][Pi[1]][Pi[2]];
        }
        else
        {
            Pi = kugelindex(Pi);
            if (Error == NO_ERRORS)
                return K[Pi[0]][Pi[1]][Pi[2]];
            else
            {
                Error = NOT_FOUND;
                return dummy;
            }
        }
    }
 
 
    template <class T> void SuperArray<T>::clear()
    {
        int anzx, anzx2;
        anzx = nges[0];
        anzx2 = nges[0] / 2;
     //   if (!iscleared)
        if (type == IN_OBJECT)
        {
            if (numObjs > 0)
            {
                for (int i = numObjs - 1; i >= 0; i--)
                {
                    if (G[i].size()>0)
                    {
                        for (int ix = n[i][0]; ix >= 0; ix--)
                        {
                            for (int iy = n[i][1]; iy >= 0; iy--)
                                G[i][ix][iy].clear();
                            G[i][ix].clear();                            
                        } // for ix    
                        G[i].clear();
                    } // if (G[i]!=NULL)
                } // for i
                n.clear(); 
                Pul.clear();
              /*  if (Obj != NULL)
                free(Obj);
                numObjs = 0;
                Obj = NULL;*/
                iscleared = true;
                
            }
        }
        else
        {
            if (K.size()>0)
            {
                for (int k = 0; k < anzx2; k++)
                {
                    /*for (int l = 0; l < ywerte[k]; l++)
                    {
                        delete[] K[anzx2 - 1 - k][ywerte[k] + l];
                        delete[] K[anzx2 + k][ywerte[k] + l];
                        delete[] K[anzx2 - 1 - k][ywerte[k] - 1 - l];
                        delete[] K[anzx2 + k][ywerte[k] - 1 - l];
                    }
                    */
                    //delete[] K[anzx2 - 1 - k];
                    //delete[] K[anzx2 + k];
                    zwerte[k].clear();                    
                }
 
                ywerte.clear();
                zwerte.clear();
 
               /* if (numObjs > 0)
                    free(Obj);*/
            }
            
            // numObjs = 0;
            iscleared = true;
        }
    }
 
    template<class T>
    inline void SuperArray<T>::reinit(double r0, INDEX_TYPE n)
    {
        reinit(r0, n, n, n);
    }
 
    template <class T>
    void SuperArray<T>::reinit(double new_r0,
        INDEX_TYPE nx,
        INDEX_TYPE ny,
        INDEX_TYPE nz
      )
    {
        this->nges = maths::Vector<INDEX_TYPE>(nx, ny, nz);
        // 2. Parameter übernehmen
        this->r0 = new_r0;
        // 1. Alte Daten freigeben
        this->arrangeGrids();
 
        
 
        // 4. Falls IN_HOST, das Kugel-Gitter (Host-Grid) neu allokieren
        if (this->type & IN_HOST)
        {
            this->allockugel();
        }
    }
 
 
    template <class T> void SuperArray<T>::copy(const SuperArray<T>& S)  // MUSS DRINGEND GE�NDERT WERDEN !
    {
        clear();
        addInc(S.Obj, S.numObjs);
        for (int i = 0; i < numObjs; i++)
            if (Obj[i]->isActive())
                for (INDEX_TYPE ix = 0; ix < n[i][0]; ix++)
                    for (INDEX_TYPE iy = 0; iy < n[i][1]; iy++)
                        for (INDEX_TYPE iz = 0; iz < n[i][2]; iz++)
                            G[i][ix][iy][iz] = S.G[i][ix][iy][iz];
        K = S.K;
        Pul = S.Pul;
        n = S.n;
        nges = S.nges;
        d = S.d;
        r0 = S.r0;
        isequal = S.isequal;
        iscleared = S.iscleared;
        pc = S.pc;
        H = S.H;
        R = S.R;
        
    }
 
    template <class T> SuperArray<T>& SuperArray<T>::operator = (const SuperArray<T>& S)
    {
        int i = 0;
        int anzx, anzx2;
 
 
        if (this == &S) return *this;
        clear();
        if (this->numObjs > 0)
        {
        /*  for (int i = 0; i < this->numObjs; i++)
                delete this->Obj[i];*/
            this->Obj.clear();
            this->Obj.shrink_to_fit();
            this->numObjs = 0;
        }
        r0 = S.r0;
        d = S.d;
        nges = S.nges;
        type = S.type;
 
        anzx = nges[0]; anzx2 = nges[0] / 2;
 
        if (type == IN_OBJECT)
        {
            isequal = true; addInc(S.Obj, S.numObjs, true);
            isequal = true;
            G = S.G;
          /*  if (S.numObjs != 0)
                do
                {
                    for (int ix = 0; ix < n[i][0]; ix++)
                        for (int iy = 0; iy < n[i][1]; iy++)
                            for (int iz = 0; iz < n[i][2]; iz++)
                                G[i][ix][iy][iz] = S.G[i][ix][iy][iz];
                    i++;
                } while (i < S.numObjs);*/
                isequal = false;
        }
        else
        {
            allockugel();
            isequal = true; addInc(S.Obj, S.numObjs, true);
            isequal = true;
            
            for (int k = 0; k < anzx2; k++)
            {
                ywerte[k] = S.ywerte[k];                            
                for (int l = 0; l < S.ywerte[k]; l++)
                {
                    zwerte[k][l] = S.zwerte[k][l];
                    for (int m = 0; m < S.zwerte[k][l]; m++)
                    {
                        K[anzx2 - 1 - k][ywerte[k] - 1 - l][m] = S.K[anzx2 - 1 - k][ywerte[k] - 1 - l][m];
                        K[anzx2 + k][ywerte[k] - 1 - l][m] = S.K[anzx2 + k][ywerte[k] - 1 - l][m];
                        K[anzx2 - 1 - k][ywerte[k] + l][m] = S.K[anzx2 - 1 - k][ywerte[k] + l][m];
                        K[anzx2 + k][ywerte[k] + l][m] = S.K[anzx2 + k][ywerte[k] + l][m];
                        K[anzx2 - 1 - k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m] = S.K[anzx2 - 1 - k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m];
                        K[anzx2 + k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m] = S.K[anzx2 + k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m];
                        K[anzx2 - 1 - k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m] = S.K[anzx2 - 1 - k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m];
                        K[anzx2 + k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m] = S.K[anzx2 + k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m];
                    }
                }
            }
        }
        return *this;
    }
 
 
 
    template <class T> void SuperArray<T>::add(const SuperArray<T>& S)
    {
        int anzx2 = nges[0] / 2;
        isequal = false;
        H = S.H;
        R = S.R;
        if (type == IN_OBJECT)
        {
            for (int i = 0; i < S.numObjs; i++)
                for (int ix = 0; ix < S.n[i][0]; ix++)
                    for (int iy = 0; iy < S.n[i][1]; iy++)
                        for (int iz = 0; iz < S.n[i][2]; iz++)
                            G[i][ix][iy][iz] += S.G[i][ix][iy][iz];
        }
        else
        {
            for (int k = 0; k < anzx2; k++)
            {
                for (int l = 0; l < ywerte[k]; l++)
                {
 
                    for (int m = 0; m < zwerte[k][l]; m++)
                    {
                        K[anzx2 - 1 - k][ywerte[k] - 1 - l][m] += S.K[anzx2 - 1 - k][ywerte[k] - 1 - l][m];
                        K[anzx2 + k][ywerte[k] - 1 - l][m] += S.K[anzx2 + k][ywerte[k] - 1 - l][m];
                        K[anzx2 - 1 - k][ywerte[k] + l][m] += S.K[anzx2 - 1 - k][ywerte[k] + l][m];
                        K[anzx2 + k][ywerte[k] + l][m] += S.K[anzx2 + k][ywerte[k] + l][m];
                        K[anzx2 - 1 - k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m] += S.K[anzx2 - 1 - k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m];
                        K[anzx2 + k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m] += S.K[anzx2 + k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m];
                        K[anzx2 - 1 - k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m] += S.K[anzx2 - 1 - k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m];
                        K[anzx2 + k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m] += S.K[anzx2 + k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m];
                    }
                }
            }
        }
    }
 
    template <class T> void SuperArray<T>::sub(const SuperArray<T>& S)
    {
        int anzx2 = nges[0] / 2;
 
        isequal = false;
        if (type == IN_OBJECT)
        {
            for (int i = 0; i < S.numObjs; i++)
                for (int ix = 0; ix < S.n[i][0]; ix++)
                    for (int iy = 0; iy < S.n[i][1]; iy++)
                        for (int iz = 0; iz < S.n[i][2]; iz++)
                            G[i][ix][iy][iz] -= S.G[i][ix][iy][iz];
        }
        else
        {
            for (int k = 0; k < anzx2; k++)
            {
                for (int l = 0; l < ywerte[k]; l++)
                {
 
                    for (int m = 0; m < zwerte[k][l]; m++)
                    {
                        K[anzx2 - 1 - k][ywerte[k] - 1 - l][m] -= S.K[anzx2 - 1 - k][ywerte[k] - 1 - l][m];
                        K[anzx2 + k][ywerte[k] - 1 - l][m] -= S.K[anzx2 + k][ywerte[k] - 1 - l][m];
                        K[anzx2 - 1 - k][ywerte[k] + l][m] -= S.K[anzx2 - 1 - k][ywerte[k] + l][m];
                        K[anzx2 + k][ywerte[k] + l][m] -= S.K[anzx2 + k][ywerte[k] + l][m];
                        K[anzx2 - 1 - k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m] -= S.K[anzx2 - 1 - k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m];
                        K[anzx2 + k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m] -= S.K[anzx2 + k][ywerte[k] - 1 - l][2 * zwerte[k][l] - 1 - m];
                        K[anzx2 - 1 - k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m] -= S.K[anzx2 - 1 - k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m];
                        K[anzx2 + k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m] -= S.K[anzx2 + k][ywerte[k] + l][2 * zwerte[k][l] - 1 - m];
                    }
                }
            }
        }
    }
 
    template <class T> std::ostream& operator << (std::ostream& os, const SuperArray<T>& S)
    {
        os << "r0=" << S.r0 << std::endl;
        os << "Ausdehnung:    nx=" << S.nges[0] << "  ny=" << S.nges[1] << "  nz=" << S.nges[2] << std::endl;
        os << S.numObjs << " Einschluesse" << std::endl;
        if (S.numObjs > 0)
            for (int i = 0; i < S.numObjs; i++)
            {
                os << "====================== Einschluss Nr. " << i << " ======================" << std::endl;
                for (int ix = 0; ix < S.n[i][0]; ix++)
                    for (int iy = 0; iy < S.n[i][1]; iy++)
                        for (int iz = 0; iz < S.n[i][2]; iz++)
                            os << abs(S.G[i][ix][iy][iz]) << std::endl;
            }
        return os;
    } 
  
    std::ostream& operator << (std::ostream& os, const SuperArray<std::vector<gridEntry> >& S);
    
 
    template <class T> void SuperArray<T>::fill(const T& x)
    {
        int anzx, anzx2;
        anzx = nges[0]; anzx2 = nges[0] / 2;
 
        if (type == IN_OBJECT)
        {
            for (int i = 0; i < numObjs; i++)
                if (Obj[i]->isActive())
                for (int ix = 0; ix < n[i][0]; ix++)
                    for (int iy = 0; iy < n[i][1]; iy++)
                        for (int iz = 0; iz < n[i][2]; iz++)
                            G[i][ix][iy][iz] = x;
        }
        else
        {
            for (int k = 0; k < anzx2; k++)
            {
                for (int l = 0; l < ywerte[k]; l++)
                {
                    for (int m = 0; m < zwerte[k][l]; m++)
                    {
                        K[anzx2 - 1 - k][ywerte[k] - 1 - l][zwerte[k][l] - 1 - m] = x;
                        K[anzx2 + k][ywerte[k] - 1 - l][zwerte[k][l] - 1 - m] = x;
                        K[anzx2 - 1 - k][ywerte[k] + l][zwerte[k][l] - 1 - m] = x;
                        K[anzx2 + k][ywerte[k] + l][zwerte[k][l] - 1 - m] = x;
                        K[anzx2 - 1 - k][ywerte[k] - 1 - l][zwerte[k][l] + m] = x;
                        K[anzx2 + k][ywerte[k] - 1 - l][zwerte[k][l] + m] = x;
                        K[anzx2 - 1 - k][ywerte[k] + l][zwerte[k][l] + m] = x;
                        K[anzx2 + k][ywerte[k] + l][zwerte[k][l] + m] = x;
                    }
                }
            }
        }
    }
 
 
    template <class T> void SuperArray<T>::makeReal()
    {
        for (int i = 0; i < numObjs; i++)
            for (int ix = 0; ix < n[i][0]; ix++)
                for (int iy = 0; iy < n[i][1]; iy++)
                    for (int iz = 0; iz < n[i][2]; iz++)
                        for (int j = 0; j < 3; j++)
                            G[i][ix][iy][iz][j] = abs(G[i][ix][iy][iz][j]);
    }
 
 
    
 
    template <class T> void SuperArray<T>::allockugel()
    {
        int anzx, anzx2;
        double dx, dy, dz, hilf;
        anzx = nges[0];
        anzx2 = nges[0] / 2;
        // cout << "allockugel" << std::endl;
      
        /*dx=d[0];
        dy=d[1];
        dz=d[2];*/
        dx = 2.0 / nges[0];
        dy = 2.0 / nges[1];
        dz = 2.0 / nges[2];
 
        zwerte.resize(anzx2);
        for (int i = 0; i < anzx2; i++)
        {
            ywerte.push_back(ceil(sqrt(1.0 - (i * dx) * (i * dx)) / dy));
            for (int j = 0; j < anzx2; j++)
            {
                hilf = 1.0 - (i * dx) * (i * dx) - (j * dy) * (j * dy);
                //      cout << "hilf:" << hilf << std::endl;
                if (hilf >= 0)
                    //        zwerte[i][j]=int(ceil(sqrt(1.0-(i*dx)*(i*dx)-(j*dy)*(j*dy))*nges[2]/2.0));
                    zwerte[i].push_back(ceil(sqrt(hilf) / dz));
                else
                    zwerte[i].push_back(1);
 
                //      cout << "zwerte[" << i << "][" << j << "]:" <<  zwerte[i][j] << std::endl;
            }
        }
 
        K.resize(anzx);
        for (int k = 0; k < anzx2; k++)
        {
            K[anzx2 - 1 - k].resize(2 * ywerte[k]);
            K[anzx2 + k].resize(2 * ywerte[k]);
            for (int l = 0; l < ywerte[k]; l++)
            {
                K[anzx2 - 1 - k][ywerte[k] + l].resize(2 * zwerte[k][l]);                
                K[anzx2 + k][ywerte[k] + l].resize(2 * zwerte[k][l]);
                K[anzx2 - 1 - k][ywerte[k] - 1 - l].resize(2 * zwerte[k][l]);
                K[anzx2 + k][ywerte[k] - 1 - l].resize(2 * zwerte[k][l]);                
            }
        }
    }
 
    template <class T> maths::Vector<INDEX_TYPE> SuperArray<T>::kugelindex(maths::Vector<INDEX_TYPE> Pi)
    {
        int anzx, jx, jy, jz, ix, iy, iz, jxh, jyh;
        maths::Vector<INDEX_TYPE> pk;
        anzx = nges[0];
        ix = Pi[0]; iy = Pi[1]; iz = Pi[2];
        Error = NO_ERRORS;
        jx = ix;
        jxh = abs(int(jx - (anzx - 1) / 2.0));
        jy = iy - (anzx / 2 - ywerte[jxh]);
        if ((jy < 0) || (jy > (2 * ywerte[jxh] - 1)))   // Punkt ausserhalb 
        {
            Error = SUPERGITTER;
            //    cout << "ix:" << ix << ", iy:" << iy << ", iz:" << iz <<  ", jx:" << jx <<
            //             ", jxh:" << jxh << ", jy:" << jy << std::endl;
            //    cout << "ywerte[" << jxh << "]:" << ywerte[jxh] << std::endl;
            return maths::Vector<INDEX_TYPE>(-1, -1, -1);
        }
        else
        {
            jyh = abs(INDEX_TYPE(jy - (2 * ywerte[jxh] - 1) / 2.0));
            jz = iz - (anzx / 2 - zwerte[jxh][jyh]);
 
            if ((jz < 0) || (jz > (2 * zwerte[jxh][jyh]) - 1)) // Punkt ausserhalb 
            {
                Error = SUPERGITTER;
                //    cout << "ix:" << ix << ", iy:" << iy << ", iz:" << iz << ", jx:" << jx << ", jy:" << jy << 
                //             ", jxh:" << jxh << ", jyh:" << jyh << ", jz:" << jz << std::endl;
                //    cout << "ywerte[" << jxh << "]:" << ywerte[jxh] << ", zwerte[" << jxh << "][" << jyh << "]:" <<
                //    zwerte[jxh][jyh] << std::endl;
                return maths::Vector<INDEX_TYPE>(-1, -1, -1);
                return maths::Vector<INDEX_TYPE>(-1, -1, -1);
            }
            else
            {
                return  maths::Vector<INDEX_TYPE>(jx, jy, jz);
            }
        }
        return pk;
    }
 
    template<class T> T SuperArray<T>::kugelwert(maths::Vector<int> Pi)
    {
        int anzx, jx, jy, jz, ix, iy, iz, jxh, jyh;
        // maths::Vector<std::complex<double> > pc;
        anzx = nges[0];
        ix = Pi[0]; iy = Pi[1]; iz = Pi[2];
        jx = ix;
        jxh = abs(int(jx - (anzx - 1) / 2.0));
        jy = iy - (anzx / 2 - ywerte[jxh]);
        pc = maths::Vector<std::complex<double> >(0.0, 0.0, 0.0);
        if ((jy < 0) || (jy > (2 * ywerte[jxh] - 1)))
        {
            std::cout << "Fehler!! Index iy falsch" << std::endl;
            return pc;
        }
        else
        {
            jyh = abs(int(jy - (2 * ywerte[jxh] - 1) / 2.0));
            jz = iz - (anzx / 2 - zwerte[jxh][jyh]);
 
            if ((jz < 0) || (jz > (2 * zwerte[jxh][jyh]) - 1))
            {
                return pc;
            }
            return K[jx][jy][jz];
        }
        return pc;
    }
 
    template<class T> T SuperArray<T>::kugelwert(int ix, int iy, int iz)
    {
        int anzx, jx, jy, jz, jxh, jyh;
        T pc;
        anzx = nges[0];
 
 
        jx = ix;
        jxh = abs(int(jx - (anzx - 1) / 2.0));
        jy = iy - (anzx / 2 - ywerte[jxh]);
        if ((jy < 0) || (jy > (2 * ywerte[jxh] - 1)))
        {
            std::cout << "Fehler!! Index iy falsch" << std::endl;
           // pc = maths::Vector<std::complex<double> >(0.0, 0.0, 0.0);
            return pc;
        }
        else
        {
            jyh = abs(int(jy - (2 * ywerte[jxh] - 1) / 2.0));
            jz = iz - (anzx / 2 - zwerte[jxh][jyh]);
            if ((jz < 0) || (jz > (2 * zwerte[jxh][jyh]) - 1))
            {
                std::cout << "Fehler!! Index iz falsch" << std::endl;
               // pc = maths::Vector<std::complex<double> >(0.0, 0.0, 0.0);
                return pc;
            }
            pc = K[jx][jy][jz];
        }
        return pc;
    }
 
    
    bool saveExPhase(SuperArray<maths::Vector<std::complex<double> > > &S, char* FName, int i = 0); 

    bool saveEyPhase(SuperArray<maths::Vector<std::complex<double> > > &S, char* FName, int i = 0);

    bool saveEzPhase(SuperArray<maths::Vector<std::complex<double> > > &S, char* FName, int i = 0);    
    bool saveExPol(SuperArray < maths::Vector < std::complex<double> > > &S, char* FName, int i = 0);
    bool saveEyPol(SuperArray < maths::Vector < std::complex<double> > > &S, char* FName, int i = 0);
    bool saveEzPol(SuperArray < maths::Vector < std::complex<double> > > &S, char* FName, int i = 0);
    bool saveabsE(SuperArray < maths::Vector < std::complex<double> > > &S, std::string FName, int i = 0);
    bool saveabsEbin(SuperArray < maths::Vector < std::complex<double> > >& S, std::string FName, int i = 0); 
    bool saveFullE(SuperArray < maths::Vector < std::complex<double> > > &S, std::string FName, int i = 0);
    double sumabs(const SuperArray<maths::Vector<std::complex<double> > >& S);
    double sumabs2(const SuperArray<maths::Vector<std::complex<double> > >& S);
    double sumabs2(const SuperArray<maths::Vector<std::complex<double> > >& S);
    double sumabs2(const SuperArray<maths::Vector<std::complex<double> > >& S, int i);
    bool save(SuperArray<GOAT::raytracing::gridEntry > S, std::string FName);
    double abs2sum(const SuperArray<maths::Vector<std::complex<double> > >& S);
bool save(SuperArray<std::vector<GOAT::raytracing::gridEntry > > S, std::string FName);      
      }
}