/************************************************************************
* Name:        Module 020 Apply Edge Rounding to Structure
* Description: The module can be used to apply rounded edges to a data array that
*              represents the height profile of a sampled interface.
*              The generated data array can be used within the optical setup for analyzing the
*              optical effect of the fabrication tolerance.
* Authors:     Christian Hellmann (LightTrans)
* Version:     1.1 (December 17th, 2014)
* License:     Some rights reserved, CC-BY 3.0
*************************************************************************/
using System;
using System.Collections.Generic;
using System.Drawing;

using VirtualLab.Programming;
using VirtualLabAPI.Core.Common;
using VirtualLabAPI.Core.FieldRepresentations;
using VirtualLabAPI.Core.Numerics;
using VirtualLabAPI.Core.Modules;

public class VLModule : IVLModule {
    public void Run() {
        //set up the internal parameters
        DataArray2D dataOfSampledInterface = null;
        double edgeRadius = 0;
        double samplingDistance = 0;
        
        #region user input
        int indexDummy = 0;
        //data matrix that shall be manipulated
        dataOfSampledInterface = (DataArray2D)(Globals.ActiveDocumentHistory.BrowseLastDocuments(DocumentFilter.DataArray2D,
                                                                                                 out indexDummy,
                                                                                                 "Select the data which shall be manipulated"));
        //check whether data array was selected
        if(dataOfSampledInterface == null){
            Globals.DataDisplay.LogError("No data field selected.");
            return;
        }
        
        //input for edge radius via GUI element
        edgeRadius = VL_GUI.AskForDouble(2000e-9,                                   //default value
                                        "Enter the Edge Radius to be applied",      //caption of the dialog
                                        "Edge Radius",                              //name of the variable to input
                                        1e-12,                                      //minimum value
                                        1e12,                                       //maximum value
                                        PhysicalProperty.Length);                   //physical property of the input parameter (length)

        if(double.IsNaN(edgeRadius)){
            Globals.DataDisplay.LogError("Canceled by user.");
            return;
        }
        //input for edge radius via GUI element
        samplingDistance = VL_GUI.AskForDouble(200e-9,                                     //default value
                                               "Enter the Sampling Distance to be used",   //caption of the dialog
                                               "Sampling Distance",                        //name of the variable to input
                                               1e-15,                                      //minimum value
                                               1e12,                                       //maximum value
                                               PhysicalProperty.Length);                   //physical property of the input parameter (length)
        if(double.IsNaN(samplingDistance)){
            Globals.DataDisplay.LogError("Canceled by user.");
            return;
        }
        #endregion
        
        #region consistency checks for user input
        //check number of subsets
        if(dataOfSampledInterface.Data.Length != 1){
            Globals.DataDisplay.LogError("The selected data array contains more or less than one subset.");
            return;
        }
        
        //check whether subsets are real-valued
        if(dataOfSampledInterface.IsComplex){
            Globals.DataDisplay.LogError("The selected data array contains a complex subset.");
            return;
        }
        
        //check whether data array is equidistance
        if(!dataOfSampledInterface.Equidistant_X ||
           !dataOfSampledInterface.Equidistant_Y){
            Globals.DataDisplay.LogError("The selected data array is not equidistant in x- and y- direction.");   
        }
        
        if(dataOfSampledInterface.InterpolationMethodForEquidistantSampling_X != dataOfSampledInterface.InterpolationMethodForEquidistantSampling_Y){
            Globals.DataDisplay.LogError("The interpolation methods in x- and y-direction do not match.");
            return;
        } 
        #endregion
        
        //extract min/max from original height values and calculate height necessary for scaling at the end
        ComplexField cfDataOriginal = dataOfSampledInterface.Data[0];
        double minHeight = 0;
        double maxHeight = 0;
        //extract min max values of the height profile
        ComplexFieldEvaluation.GetMinMaxValues(cfDataOriginal,
                                               FieldQuantity.RealPart,
                                               out minHeight,
                                               out maxHeight);
        
        //evaluate sampling parameters for convolution
        SamplingParameters sParaConvolution = new SamplingParameters();
        
        VectorD diameterDataArray = new VectorD(dataOfSampledInterface.GetCoordinateRange_X(dataOfSampledInterface.InterpolationMethodForEquidistantSampling_X == InterpolationMethod.Nearest).SignedValueRange,
                                                dataOfSampledInterface.GetCoordinateRange_Y(dataOfSampledInterface.InterpolationMethodForEquidistantSampling_Y == InterpolationMethod.Nearest).SignedValueRange);
        //calculate number of sampling points
        sParaConvolution.SamplingPoints = (Vector)(diameterDataArray / new VectorD(samplingDistance, samplingDistance));
        
        //ensure that the result has always odd number of data points
        if(sParaConvolution.SamplingPoints.X % 2 == 0){
            sParaConvolution.SamplingPoints += new Vector(1,0);
        }
        if(sParaConvolution.SamplingPoints.Y % 2 == 0){
            sParaConvolution.SamplingPoints += new Vector(0,1);
        }
        
        //calculate samping distance
        sParaConvolution.SamplingDistance = diameterDataArray / sParaConvolution.SamplingPoints ;
        
        //generate gaussian beam with correct size and sampling
        double tempWavelength = 500e-9; //we need to specify a wavelength => only temporary
        //calculate waist from radius (user input)
        VectorD waist = new VectorD(edgeRadius / 2.0,
                                    edgeRadius / 2.0); 
        
        //call VirtualLab method to generate a gaussian function 
        ComplexAmplitude caGaussian = WaveFieldGenerator.GaussianWave(sParaConvolution,        //sampling parameters for gaussian
                                                                      waist,                   //the waist (calculated from radius + temp wavelength
                                                                      new VectorD(0, 0),       //longitudinal position
                                                                      new VectorD(0, 0),       //lateral position
                                                                      new Vector(0,0),         //order
                                                                      tempWavelength,          //wavelength (only used temporary)
                                                                      true,                    //is hermite?
                                                                      PrecisionMode.Double,    //precision mode of the field which shall be generated
                                                                      false,                   //is complex?
                                                                      true);                   //propagate in positiv z-direction (temporary)
        
        ComplexField cfGaussian = caGaussian.Field;
                
        #region interpolate data of sampled interface to correct sampling   
        //extract sampling parameters from original data array
        VectorD samplingDistanceOriginalField = new VectorD(dataOfSampledInterface.SamplingDistance_X,
                                                            dataOfSampledInterface.SamplingDistance_Y);
        Vector numberSamplingPointsOriginalField = new Vector(dataOfSampledInterface.NoOfDataPoints_X,
                                                              dataOfSampledInterface.NoOfDataPoints_Y);  
           
        //calculate zero point
        VectorD centerDataPoint = new VectorD(dataOfSampledInterface.CoordinateOfFirstDataPoint_X + (dataOfSampledInterface.SamplingDistance_X * dataOfSampledInterface.NoOfDataPoints_X / 2.0),
                                              dataOfSampledInterface.CoordinateOfFirstDataPoint_Y + (dataOfSampledInterface.SamplingDistance_Y * dataOfSampledInterface.NoOfDataPoints_Y / 2.0));

        ComplexField[] fields = new ComplexField[dataOfSampledInterface.DimensionalityOfData];
        //calculate sampling parameters for pixel interpolation
        SamplingParameters sParaPixel = new SamplingParameters(sParaConvolution.SamplingPoints,
                                                               sParaConvolution.SamplingDistance / samplingDistanceOriginalField);
        //calculate new center point in pixel coordinates
        VectorD centerPointPixel = (numberSamplingPointsOriginalField / 2) -
                                   ((sParaPixel.SamplingPoints / 2) * sParaPixel.SamplingDistance);

        //special handling for even number of data points => shift by -0.5 Pixel
        if((numberSamplingPointsOriginalField.X % 2) == 0){
            centerPointPixel.X -= 0.5;
        }
        if((numberSamplingPointsOriginalField.Y % 2) == 0){
            centerPointPixel.Y -= 0.5;
        }
                
        //perform interpolation to data fields
        for (int runDataSubSets = 0; runDataSubSets < dataOfSampledInterface.DimensionalityOfData; runDataSubSets++) {
            fields[runDataSubSets] = ComplexFieldInterpolation.Interpolation(dataOfSampledInterface.Data[runDataSubSets],
                                                                             sParaPixel,
                                                                             centerPointPixel,
                                                                             dataOfSampledInterface.InterpolationMethodForEquidistantSampling_X);
        }

        //perform interpolation of Coordinates
        double samplingDistanceX = sParaConvolution.SamplingDistanceX;
        double samplingDistanceY = sParaConvolution.SamplingDistanceY;
        //calculate position of first pixel ==> it is always ensured that the data matrix is centered around (0,0)
        double coordinateOfFirstDataPointX = ((-1) * (sParaConvolution.Diameter.X / 2.0)) +
                                             (sParaConvolution.SamplingDistance.X / 2.0);
        double coordinateOfFirstDataPointY = ((-1) * (sParaConvolution.Diameter.Y / 2.0)) +
                                             (sParaConvolution.SamplingDistance.Y / 2.0);
        
        //construct new data array which contains the interpolated data
        DataArray2D interpolatedFabricationData = new DataArray2D(new ComplexFieldArray(fields),
                                                                  dataOfSampledInterface.PhysicalPropertiesOfDataEntries,
                                                                  dataOfSampledInterface.CommentsOfDataEntries,
                                                                  samplingDistanceX,
                                                                  coordinateOfFirstDataPointX,
                                                                  dataOfSampledInterface.PhysicalPropertyOf_X_Coordinates,
                                                                  dataOfSampledInterface.CommentOfCoordinates_X,
                                                                  samplingDistanceY,
                                                                  coordinateOfFirstDataPointY,
                                                                  dataOfSampledInterface.PhysicalPropertyOf_Y_Coordinates,
                                                                  dataOfSampledInterface.CommentOfCoordinates_Y);
        #endregion
        //extract field from interpolated data array
        ComplexField cfFabricationData = interpolatedFabricationData.Data[0];
        
        //perform convolution
        VirtualLabAPI.Core.Numerics.Transformations.Fourier(cfGaussian);
        VirtualLabAPI.Core.Numerics.Transformations.Fourier(cfFabricationData);
        cfFabricationData.MultEqual(cfGaussian);
        VirtualLabAPI.Core.Numerics.Transformations.InverseFourier(cfFabricationData);
        
        //ensure that no complex - valued fields are generated
        cfFabricationData.IsComplex = false;
        //lift positive + normalize + correct scaling
        cfFabricationData.LiftPositive();
        cfFabricationData.Normalize();        
        cfFabricationData.MultEqual(maxHeight - minHeight);
        interpolatedFabricationData.Data[0] = cfFabricationData;
        
        //set correct interpolation method for data array
        interpolatedFabricationData.InterpolationMethodForEquidistantSampling_X = InterpolationMethod.Cubic8P;
        interpolatedFabricationData.InterpolationMethodForEquidistantSampling_Y = InterpolationMethod.Cubic8P;
        
        //output of modified data array
        VL_GUI.ShowDocument(interpolatedFabricationData, "Fabrication Data with Edge Tolerances");
    }
}
