ibug aus calcArg raus

program2D/source/extemp/calcArg.cc

@@ -74,8 +74,16 @@ CalcContinuousArg::CalcContinuousArg(Blockgrid2D& blockgrid_) {
     
     orderedPoints = new std::list<PointForPhase>[numberBlocks];
         
-    enum woSmaller { nichts, Norden, Sueden, Osten, Westen };
+    enum woSmaller { nichts, Norden, Sueden, Osten, Westen, NOst, SOst, NWest, SWest };
         
+    idOrigin = 0;
+    iOrigin  = 0;
+    jOrigin  = 0;
+    NxOrigin = 0;
+    
+    bool nothingFound    = true;
+    double smallesRadius = 0.0;
+    
     for(int id=0;id<numberBlocks;++id) {
         int Nx = blockgrid->Give_Nx_rectangle(id);
         int Ny = blockgrid->Give_Ny_rectangle(id);
@@ -155,13 +163,70 @@ CalcContinuousArg::CalcContinuousArg(Blockgrid2D& blockgrid_) {
 		  woNext = Norden;
 	       }
 	    }
-
+	        
+	    //Nord-Westen
+            if(i>0 && j<Ny) { 
+	       double r = norm.template Give_data<rectangleEl>(id,i-1,j+1,Nx);
+               if(r < radiusNext) {
+		  radiusNext = r;
+		  woNext = NWest;
+	       }
+	    }
+	    //Nord-Osten
+            if(i<Nx && j<Ny) { 
+	       double r = norm.template Give_data<rectangleEl>(id,i+1,j+1,Nx);
+               if(r < radiusNext) {
+		  radiusNext = r;
+		  woNext = NOst;
+	       }
+	    }
+	    //Sued-Westen
+            if(i>0 && j>0) { 
+	       double r = norm.template Give_data<rectangleEl>(id,i-1,j-1,Nx);
+               if(r < radiusNext) {
+		  radiusNext = r;
+		  woNext = SWest;
+	       }
+	    }
+	    //Sued-Osten
+            if(i<Nx && j>0) { 
+	       double r = norm.template Give_data<rectangleEl>(id,i+1,j-1,Nx);
+               if(r < radiusNext) {
+		  radiusNext = r;
+		  woNext = SOst;
+	       }
+	    }
+	    
+	    if(woNext == nichts) {
+	       if(nothingFound) {
+		  nothingFound = false;
+		  idOrigin = id;
+                  iOrigin  = i;
+                  jOrigin  = j;
+		  NxOrigin = Nx;
+                  smallesRadius = radiusM;
+	       }
+	       else {
+		  if(smallesRadius > radiusM) {
+		     idOrigin = id;
+                     iOrigin  = i;
+                     jOrigin  = j;
+  		     NxOrigin = Nx;		     
+                     smallesRadius = radiusM;		    
+		  }
+	       }
+	    }
 
 
 	    if(woNext == Westen) orderedPoints[id].push_back(PointForPhase(i,j,radiusM,i-1,j));
 	    if(woNext == Osten)  orderedPoints[id].push_back(PointForPhase(i,j,radiusM,i+1,j));
 	    if(woNext == Sueden) orderedPoints[id].push_back(PointForPhase(i,j,radiusM,i,j-1));
-	    if(woNext == Norden) orderedPoints[id].push_back(PointForPhase(i,j,radiusM,i,j+1));   
+	    if(woNext == Norden) orderedPoints[id].push_back(PointForPhase(i,j,radiusM,i,j+1));
+	    
+	    if(woNext == NWest) orderedPoints[id].push_back(PointForPhase(i,j,radiusM,i-1,j+1));
+	    if(woNext == NOst)  orderedPoints[id].push_back(PointForPhase(i,j,radiusM,i+1,j+1));	    
+	    if(woNext == SWest) orderedPoints[id].push_back(PointForPhase(i,j,radiusM,i-1,j-1));
+	    if(woNext == SOst)  orderedPoints[id].push_back(PointForPhase(i,j,radiusM,i+1,j-1));	    	    
 	}
 	orderedPoints[id].sort(compare_Radius);    
 	
@@ -181,10 +246,17 @@ double phaseComp(std::complex<double> value) {
   return atan2(value.imag(),value.real());
 }
 
-void CalcContinuousArg::calcArg(Variable2D<double>& phase, Variable2D<std::complex<double> >& Efield) {
+void CalcContinuousArg::calcArg(Variable2D<double>& phase, Variable2D<std::complex<double> >& Efield, bool zeroAtOrigin) {  
+   double phaseOrigin = 0.0;
+   
+   if(zeroAtOrigin) {
+      std::complex<double> EField_origin = Efield.template Give_data<rectangleEl>(idOrigin,iOrigin,jOrigin,NxOrigin);   
+      phaseOrigin = atan2(EField_origin.imag(),EField_origin.real());        
+   }
+   
    Function2d1<double, std::complex<double> >  Phase(phaseComp);
-   phase = Phase(Efield);           
-  
+   phase = Phase(Efield) - phaseOrigin;      
+   
    for(std::list<BlockForPhase>::iterator it=orderedBlocks.begin(); it != orderedBlocks.end(); ++it) {
        int id = it->getId();
        phase.Update<rectangleEl>(id);       
@@ -195,18 +267,25 @@ void CalcContinuousArg::calcArg(Variable2D<double>& phase, Variable2D<std::compl
 	   int i = itP->getI();
 	   int j = itP->getJ();
 	   
-	   std::complex<double> EField_M = Efield.template Give_data<rectangleEl>(id,i,j,Nx);
-            
-           double phaseM = atan2(EField_M.imag(),EField_M.real());
-	   
+	   std::complex<double> EField_M = Efield.template Give_data<rectangleEl>(id,i,j,Nx); // muss man das echt nochmal holen???
+           double phaseM = atan2(EField_M.imag(),EField_M.real()) - phaseOrigin;
+	   	   
 	   int iSmall = itP->getIsmall();
 	   int jSmall = itP->getJsmall();
            
            double phaseNext = phase.template Give_data<rectangleEl>(id,iSmall,jSmall,Nx);
  	   
+
+	   int numPi = (phaseM-phaseNext)/(2.0*M_PI);
+	   phaseM = phaseM - 2.0*M_PI * numPi;
+
 	   if(fabs(phaseM-phaseNext)>M_PI) {
-	      if(phaseM < phaseNext) while(phaseM < phaseNext && fabs(phaseM-phaseNext)>M_PI) phaseM = phaseM + 2.0*M_PI;
-	      else                   while(phaseM > phaseNext && fabs(phaseM-phaseNext)>M_PI) phaseM = phaseM - 2.0*M_PI;
+	      if(phaseM < phaseNext) {
+		while(phaseM < phaseNext && fabs(phaseM-phaseNext)>M_PI) phaseM = phaseM + 2.0*M_PI;
+	      }
+	      else {
+		while(phaseM > phaseNext && fabs(phaseM-phaseNext)>M_PI) phaseM = phaseM - 2.0*M_PI;
+	      }
               phase.setValueAtRectangularPoint(phaseM,id,i,j,Nx);	      
 	   }
        }

program2D/source/extemp/calcArg.h

@@ -63,10 +63,13 @@ class BlockForPhase {
 
 /** \addtogroup ExpressionTemplates2D **/
 /* @{ */   
+/**
+ * Calculates the continuous Phase of a complex 2D Beam
+ **/
 class CalcContinuousArg {
   public:  
     CalcContinuousArg(Blockgrid2D& blockgrid_);
-    void calcArg(Variable2D<double>& phase, Variable2D<std::complex<double> >& Efield);    
+    void calcArg(Variable2D<double>& phase, Variable2D<std::complex<double> >& Efield, bool zeroAtOrigin = false);    
     
   private:
     Blockgrid2D*        blockgrid;
@@ -75,7 +78,10 @@ class CalcContinuousArg {
     std::list<BlockForPhase> orderedBlocks; 
     int numberBlocks;
    
-    
+    int idOrigin;
+    int iOrigin;
+    int jOrigin;
+    int NxOrigin;
 };
 /* @} */