//******************************************************************************
//** SCATMECH: Polarized Light Scattering C++ Class Library
//** 
//** File: grating.cpp
//**
//** Thomas A. Germer
//** Optical Technology Division, National Institute of Standards and Technology
//** 100 Bureau Dr. Stop 8443; Gaithersburg, MD 20899-8443
//** Phone: (301) 975-2876; FAX: (301) 975-6991
//** Email: thomas.germer@nist.gov
//**
//** Version: 6.01 (May 2008)
//**
//******************************************************************************
#include "grating.h"
#include <fstream>
#include <algorithm>
#include <valarray>
#include <list>
#include <limits>
#include "scateval.h"
 
using namespace std;

namespace SCATMECH {

    static inline double mymax(double x,double y) {return (x>y)?x:y;}

    void Register(const Grating* x)
    {
        static bool regd=false;
        if (!regd) {
            regd=true;

            Register_Model(Grating);

            // Future releases...
            Register_Model(Generic_Grating); 
            Register_Model(Single_Line_Grating);
            Register_Model(Sinusoidal_Grating);
            Register_Model(Triangular_Grating);
            Register_Model(Corner_Rounded_Grating);
            Register_Model(Null_Grating);
        }
    }

    Fourier_Component_Calculator::
    Fourier_Component_Calculator(double _period,int _order,int _recip,double _boundary)
            : period(_period), order(_order), sum(0.), begin(true), recip(_recip), boundary(_boundary){}
    
    void Fourier_Component_Calculator::enter(double x,COMPLEX y)
    {
        if (recip) y = 1./y;

        if (begin) {
            begin = false;
        } else {
            double x0 = last_x-boundary/2.;
            double x1 = last_x+boundary/2.;
            double x2 = x-boundary/2.;
            COMPLEX I(0,1);
            double n=order;
            COMPLEX y1=y;
            COMPLEX y0=last_y;
            if (order!=0) {
                sum += COMPLEX(0.,0.5)*
                    (exp(COMPLEX(0.,2*order*pi*x1/period)) - 
                     exp(COMPLEX(0.,2*order*pi*x2/period)))*y/(order*pi);
            
                if (boundary!=0.) {
                    double period2=sqr(period);
                    double period3=period*period2;

                    sum += (
                            COMPLEX(0,0.25)*(
                                exp((2.*I*n*pi*x0)/period)*(
                                    2.*pow(n*pi*(x0-x1),3)*y0 + 
                                    COMPLEX(0,3)*pow(period,3)*(y0 - y1) + 
                                    3.*n*pow(period,2)*pi*(x0 - x1)*(y0 - y1)
                                                           ) + 
                                exp((2.*I*n*pi*x1)/period)*(
                                    COMPLEX(0,-3.)*pow(period,3)*(y0 - y1) + 
                                    3.*n*pow(period,2)*pi*(x0 - x1)*(y0 - y1) 
                                    - 2.*pow(n*pi*(x0 - x1),3)*y1
                                                           )
                                            )
                            )/(pow(n*pi,4)*pow(x0 - x1,3));
                }
            } else {
                sum += y*(x2 - x1)/period;
                if (boundary!=0.) {
                    sum += -((x0 - x1)*(y0 + y1))/(2.*period);
                }
            }
        }
        last_x = x;
        last_y = y;
    }

    COMPLEX Grating::epsilon(const dielectric_function& e) 
    {
        return (COMPLEX)e.epsilon(lambda);
    }

    void Single_Line_Grating::setup() 
    {
        Grating::setup();
        eline = epsilon(material);
        espace = epsilon(space);
    }

    COMPLEX Single_Line_Grating::fourier(int order,int level,int recip) 
    {    
        SETUP();
    
        if (level<0||level>=nlevels) 
            error("Invalid level in fourier()");
    
        double h = (0.5+level)/(double)nlevels;
        double x1 = -topwidth/2.-(bottomwidth-topwidth)*h/2.-offset*h;
        double x2 = topwidth/2.+(bottomwidth-topwidth)*h/2.-offset*h;
    
        Fourier_Component_Calculator z(period,order,recip,boundary);
        z.enter(-period/2,espace);
        z.enter(x1,espace);
        z.enter(x2,eline);
        z.enter(period/2.,espace);
        return z.result();
    }

    void Corner_Rounded_Grating::setup()
    {
        Grating::setup();
        const char no_straight_section[]=
            "No straight section of sidewall";
        const char no_flat_section_top[]=
            "No flat section on top of line";
        const char no_flat_section_bottom[]=
            "No flat section on bottom of trench";

        eps_line=epsilon(material);

        x.resize(nlevels);
        y.resize(nlevels);

        double theta=sidewall*deg;
        double cottheta = sidewall!=90. ? 1./tan(theta) : 0.;
        double sintheta = sin(theta);
        double costheta = cos(theta);
        double r1=radiust;
        double r2=radiusb;

        if (sidewall<=90) {
            // See TAG Notes 3 April 2006...
            double a2=r1*sintheta;
            double a4=r2*sintheta;
            double a5=r2*(1-costheta)*cottheta;
            double a3=(height-(r1+r2)*(1.-costheta))*cottheta;
            double a1=width/2-a3-a2-a5;
            double a6=a2+a3+a4;

            if (a3<0) error(no_straight_section);
            if (a1<0) error(no_flat_section_top);
            if (a6+a1>period/2.) error(no_flat_section_bottom);

            for (int i=0;i<nlevels;++i) {
                double xx = (i+0.5)/(double)(nlevels)*a6;
                double yy;
                if (xx<a2) {
                    yy = height-r1+sqrt(r1*r1-xx*xx);
                } else if (xx<a2+a3) {
                    double xxx = a6-a4+a5-xx;
                    yy = xxx*tan(theta);
                } else {
                    double xxx = a6-xx;
                    yy = r2-sqrt(r2*r2-xxx*xxx); 
                }
                x[i]=a1+xx;
                y[i]=yy;
            }
        } else {
            // See TAG Notes 25 January 2007...
            double b1 = width/2+r2*((costheta-1)*cottheta-1+sintheta);
            double b2 = width/2+r2*(costheta-1)*cottheta;
            double b3 = width/2+r2*((costheta-1)*cottheta+sintheta);
            double w2 = width-2*height*cottheta;
            double c3 = w2/2-r1*((costheta-1)*cottheta-1+sintheta);
            double c2 = w2/2-r1*(costheta-1)*cottheta;
            double c1 = w2/2-r1*((costheta-1)*cottheta+sintheta);

            if (c2<b2) error(no_straight_section);
            if (c1<0) error(no_flat_section_top);
            if (b3>period/2.) error(no_flat_section_bottom);

            double tspan=c3-c1;
            double bspan=b3-b1;
            double sspan=c3-b1;
            double span=tspan+bspan+sspan;

            for (int i=0;i<nlevels;++i) {
                double xx = (i+0.5)/(double)(nlevels)*span;
                double yy;
                if (xx<tspan) {
                    double x=xx;
                    yy = height-r1+sqrt(r1*r1-x*x);
                    xx = c1+x;
                } else if (xx<tspan+c3-c2) {
                    double x = xx-tspan;
                    yy = height-r1-sqrt(sqr(r1)-sqr(r1-x));
                    xx = c3-x;
                } else if (xx<tspan+c3-b2) {
                    double x=xx-tspan-c3+c2;
                    double y1 = r2*(1-costheta);
                    double y2 = r1*(1-costheta);
                    yy = height-y2-x/(c2-b2)*(height-y1-y2);
                    xx = c2-x;
                } else if (xx<tspan+c3-b1) {
                    double x=xx-tspan-c3+b2;
                    x = b2-b1-x;
                    yy = r2+sqrt(sqr(r2)-sqr(r2-x));
                    xx = b1+x;
                } else if (xx<tspan+c3-b1+b3-b1) {
                    double x=xx-tspan-c3+b1;
                    yy = r2-sqrt(sqr(r2)-sqr(r2-x));
                    xx = b1+x;
                }
                x[i]=xx;
                y[i]=yy;
            }
        }
    }

    COMPLEX Corner_Rounded_Grating::fourier(int order,int level,int recip)
    {
        SETUP();

        if (level>=nlevels) return 0.;

        COMPLEX result = 0; 
        Fourier_Component_Calculator z(period,order,recip,boundary);
        z.enter(-period/2.,vacuum);
        z.enter(-x[level],vacuum);
        z.enter(x[level],eps_line);
        z.enter(period/2.,vacuum);
        return z.result();
    }

    COMPLEX Sinusoidal_Grating::fourier(int i,int level,int recip)
    {
        SETUP();
        int _nlevels = base!=0 ? nlevels-1 : nlevels;

        if (level==_nlevels) {
            COMPLEX e=epsilon(material);
            Fourier_Component_Calculator z(period,i,recip,Grating::boundary);
            z.enter(-period/2,e);
            z.enter(period/2,e);
            return z.result();
        } else {
            double xh = period/2.*(level+0.5)/_nlevels;
            COMPLEX vacuum(1,0);
            COMPLEX e=epsilon(material);
            Fourier_Component_Calculator z(period,i,recip,Grating::boundary);
            z.enter(-period/2,vacuum);
            z.enter(-xh,vacuum);
            z.enter(xh,e);
            z.enter(period/2,vacuum);
            return z.result();
        }
    }

    COMPLEX Triangular_Grating::fourier(int i,int level,int recip)
    {
        SETUP();

        COMPLEX vacuum(1,0);
        COMPLEX m=epsilon(material);
        double apex1=period*aspect;
        double apex2=period*(1.-aspect);
        double x1 = -period/2.+apex1*(nlevels-level-0.5)/nlevels;
        double x2 = period/2.-apex2*(nlevels-level-0.5)/nlevels;
    
        Fourier_Component_Calculator z(period,i,recip,Grating::boundary);
        z.enter(-period/2.,vacuum);
        z.enter(x1,vacuum);
        z.enter(x2,m);
        z.enter(period/2.,vacuum);
        return z.result();
    }


    namespace {
        struct segment {
            COMPLEX vertex1,vertex2;
            std::string mat1,mat2;
        };
        struct bounds {
            double x; 
            std::string mat1,mat2;

            bool operator>(const bounds& b) {
                return x>b.x;
            }
            bool operator<(const bounds& b) {
                return x<b.x;
            }

        };

        string read_pair(istream& is)
        {
            is >> ws;

            if (is.peek()!='(') {
                string result;
                is >> result;
                return result;
            }
        

            // Get the contents of the parentheses and place them in the string contents...
            string result;
            result += is.get();
            int level=1;
            while (level>0) {
                if (is.eof()) throw SCATMECH_exception("Mismatched parentheses");
                char next = is.get();
                if (next==')') --level;
                if (next=='(') ++level;
                
                if (level>=0) result += next;        
            }
            return result;
        }   

    }

    void Generic_Grating::setup()
    {
        Grating::setup();

        string fname = find_file(filename);
        ifstream_with_comments file(fname.c_str());
        if (!file) error("Cannot open file");

        int i,j;

        eps.clear();
        position.clear();
        thickness.clear();
        material.clear();
        vars.clear();
        vnames.clear();

        string stemp;

        //*****************************************************************
        //*
        //* First block: read the parameters ...
        //*
        //*****************************************************************

        // Read and check the section heading...
        file >> stemp;
        if (stemp!="PARAMETERS") error("Expected PARAMETERS label");

        // Evaluate the values passed through the parameter pstring...
        Evaluator pe(pstring);

        vars["period"] = period;
        vars["lambda"] = lambda;

        i=0;
        while (1) {
            file >> stemp;
            if (stemp=="END") break;
            if (file.fail()) error("Error reading file in PARAMETERS section");

            varsmap::iterator a = vars.find(stemp);
            if (a!=vars.end()) error("Duplicate parameter label: " + stemp);
            if (i>pe.NResult()) 
                error("Number of values in pstring (" + to_string(pe.NResult()) + 
                      ") less than number of parameters (>" + to_string(i) + ")");
            double v = pe.Result(i);
            vars.insert(a,varsmap::value_type(stemp,v));
            i++;
        }
        if (i!=pe.NResult()) 
            error("Number of values in pstring (" + to_string(pe.NResult()) + 
                  ") greater than number of parameters (" + to_string(i) + ")");


        //*****************************************************************
        //*
        //* Read in working variables...
        //*
        //*****************************************************************
        file >> stemp;
        if (stemp=="WORKING") {
            while (1) {
                file >> stemp;
                if (stemp=="END") break;
        
                string expression = file.getquoted();

                if (file.fail()) error("Error reading file in WORKING section");

                double value = (double)(Evaluator(expression,vars));

                vars[stemp]=value;
            }
            file >> stemp;
        }

        //*****************************************************************
        //*
        //* Second block: read in materials...
        //*
        //*****************************************************************

        if (stemp!="MATERIALS") error("Expected MATERIALS label");

        eps["medium_i"] = epsilon(medium_i);
        eps["medium_t"] = epsilon(medium_t);

        while (1) {
            // Read material index...
            file >> stemp;
            if (stemp=="END") break;
    
            // Read material dielectric function...
            string smaterial = read_pair(file);

            if (file.fail()) error("Error reading file in MATERIALS section");

            string smaterial2;
            try {
                smaterial2 = Evaluator(smaterial,vars).ResultString();
            } catch (SCATMECH_exception&) {
                smaterial2 = smaterial;
            }
            dielectric_function df(smaterial2);

            epsmap::iterator a = eps.find(stemp);
            if (a!=eps.end()) error("Duplicate material index: " + stemp);;
            eps.insert(a,epsmap::value_type(stemp,epsilon(df)));
        } 

        //*****************************************************************
        //*
        //* Third block: read vertices...
        //*
        //*****************************************************************

        file >> stemp;
        if (stemp!="VERTICES") error("Expected MATERIALS label");

        typedef map<string,COMPLEX> verticesmap;
        verticesmap vertices;

        while (1) {
            file >> stemp;
            if (stemp=="END") break;
            COMPLEX vertex = get_complex_value(file);
            if (file.fail()) error("Error reading file in VERTICES section");

            verticesmap::iterator a = vertices.find(stemp);
    
            if (a!=vertices.end()) error("Duplicate vertex label: " + stemp);

            vertices.insert(a,verticesmap::value_type(stemp,vertex));
        }

        //*****************************************************************
        //*
        //* Fourth block: read boundaries...
        //*
        //*****************************************************************
        file >> stemp;
        if (stemp!="BOUNDARIES") error("Expected BOUNDARIES label");

        message("Boundaries:\n");

        vector<segment> boundaries;
        while (!file.eof()) {
            string vertex1,vertex2;
            string imat1,imat2;

            // Read first vertex...
            file >> vertex1;
            if (vertex1=="END") break;
            // If end of file, then end of block...
            if (!file.eof()) {
                if (file.fail()) error("Error reading file at vertex #1");
                verticesmap::iterator v1 = vertices.find(vertex1);
                if (v1==vertices.end()) error("Unknown vertex label: " + vertex1);
            
                // Read second vertex...
                file >> vertex2;
                if (file.fail()) error("Error reading file at vertex #2");
                verticesmap::iterator v2 = vertices.find(vertex2);
                if (v2==vertices.end()) error("Unknown vertex label: " + vertex2);
            
                // Read left material index...
                file >> imat1;
                if (file.fail()) error("Error reading file at material #1");
                epsmap::iterator m1 = eps.find(imat1);
                if (m1==eps.end()) error ("Unknown material label: " + imat1);
            
                // Read right material index...
                file >> imat2;
                if (file.fail()) error("Error reading file at material #2");
                epsmap::iterator m2 = eps.find(imat2);
                if (m2==eps.end()) error ("Unknown material: " + imat2);
            
                // Load into a segment...
                segment seg;
                seg.vertex1=v1->second;
                seg.vertex2=v2->second;
                seg.mat1=imat1;
                seg.mat2=imat2;

                message(imat1 + tab + imat2 + tab + format("%g",real(seg.vertex1)) + tab + format("%g",imag(seg.vertex1)) + '\n' +
                        imat1 + tab + imat2 + tab + format("%g",real(seg.vertex2)) + tab + format("%g",imag(seg.vertex2)) + "\n\n");

                // Store in list of boundaries...
                boundaries.push_back(seg);
            }
        } 

        //*****************************************************************
        //* 
        //* Get a list of all vertical and horizontal positions...
        //*
        //*****************************************************************
        list<double> X,Y;
        for (i=0;i<boundaries.size();++i) {
            X.push_back(real(boundaries[i].vertex1));
            X.push_back(real(boundaries[i].vertex2));
            Y.push_back(imag(boundaries[i].vertex1));
            Y.push_back(imag(boundaries[i].vertex2));
        }

        // Sort them and throw away duplicates 
        X.sort();
        X.unique();
        Y.sort();
        Y.unique();

        // Total period of the grating...
        period = X.back()-X.front();

        // Total height of grating...
        double height = Y.back()-Y.front();

        double firstx = X.front();
        double prevx = firstx;
        double firsty = Y.front();
        double prevy = firsty;

        list<double>::iterator p;

        vector<double> td;
        double total_td=0;

        //*****************************************************************
        //*
        //* Check for consistency along slices in the vertical direction
        //*
        //*****************************************************************
        list<bounds> Yslice;
        for (p=++(X.begin());p!=X.end();++p) {
            // Get the horizontal position between the two x-positions
            double width = (*p-prevx);
            double x0 = prevx+width/2.;
            prevx = *p;

            Yslice.clear();
            for (j=0;j<boundaries.size();++j) {
                double x1=real(boundaries[j].vertex1);
                double y1=imag(boundaries[j].vertex1);
                double x2=real(boundaries[j].vertex2);
                double y2=imag(boundaries[j].vertex2);
                
                // If segment crosses horizontal position...
                if ((x1>=x0 && x2<x0) || (x1<=x0 && x2>x0)) {
                    bounds b;
                    // Get vertical position of crossing...
                    b.x = (x0*y1 - x0*y2 + y2*x1 - y1*x2)/(x1 - x2);

                    // Store positions and material indices...
                    if (x2>x1) {
                        b.mat1 = boundaries[j].mat2;
                        b.mat2 = boundaries[j].mat1;
                    } else {
                        b.mat1 = boundaries[j].mat1;
                        b.mat2 = boundaries[j].mat2;
                    }
                    Yslice.push_back(b);
                }
            }
            Yslice.sort();
            // Store positions and materials in microlayer
            string imat="medium_t";
            for (list<bounds>::iterator q=Yslice.begin();q!=Yslice.end();++q) {
                if (q->mat1!=imat) 
                    error("Inconsistent material index: " + q->mat1 + "!=" + imat);
                imat = q->mat2;
            }
            if (imat!="medium_i") 
                error("Inconsistent material index: " + string("medium_i") + "!=" + imat);
        }

        //*****************************************************************
        //*
        //* Find maximum traversed distance for each sublayer...
        //*
        //* Note: Sublayers refer to regions between different given 
        //*       y-values, which are further divided into microlayers.
        //*
        //*****************************************************************
        for (p=++(Y.begin());p!=Y.end();++p) {
            // Get a vertical position in the middle of the sublayer...
            double subheight=(*p-prevy);
            double h=prevy+subheight/2.;

            double max=0;
            for (j=0;j<boundaries.size();++j) {
                double x1=real(boundaries[j].vertex1);
                double y1=imag(boundaries[j].vertex1);
                double x2=real(boundaries[j].vertex2);
                double y2=imag(boundaries[j].vertex2);

                // If segment crosses vertical position...
                if ((y1>=h && y2<h) || (y1<=h && y2>h)) {
                    // Calculate the traversed distance...
                    double tdist = fabs((x2-x1)/(y2-y1)*subheight);
                    // Check to see if the materials on both sides are the same...
                    if (boundaries[j].mat1==boundaries[j].mat2) tdist=0;
                    // If this is a bigger distance, save it...
                    if (max<tdist) max=tdist;
                }
            }
            // Store maximum traversed distance...
            td.push_back(max);
            total_td += max;
            prevy=*p;
        }

        //*****************************************************************
        //*
        //* Calculate and store boundaries for each microlayer...
        //*
        //*****************************************************************
        int k=0;
        int kk=0;
        prevy=firsty;
        for (p=++(Y.begin());p!=Y.end();++p) {
            double subheight=(*p-prevy);
            // Scale number of microlayers in sublayer by maximum traversed distance....
            int mlayers = (total_td==0) ? 1 : (int)(td[kk]/total_td*nlayers+0.5);
            // Make sure there is at least one microlayer...
            if (mlayers<1) mlayers = 1;
            for (i=0;i<mlayers;++i) {
                // Calculate vertical position in middle of microlayer...
                double h=prevy+subheight*(i+0.5)/mlayers;
                // Start new microlayer...
                position.push_back(vector<double>());
                material.push_back(vector<string>());
                thickness.push_back(subheight/mlayers);

                // Get locations of boundaries...
                list<bounds> X;
                for (j=0;j<boundaries.size();++j) {
                    double x1=real(boundaries[j].vertex1);
                    double y1=imag(boundaries[j].vertex1);
                    double x2=real(boundaries[j].vertex2);
                    double y2=imag(boundaries[j].vertex2);

                    // If segment crosses vertical position...
                    if ((y1>=h && y2<h) || (y1<=h && y2>h)) {
                        bounds b;
                        // Get horizontal position of crossing...
                        b.x = (h*x1 - h*x2 + x2*y1 - x1*y2)/(y1 - y2);
                        // Store positions and material indices...
                        if (y2>y1) {
                            b.mat1 = boundaries[j].mat1;
                            b.mat2 = boundaries[j].mat2;
                        } else {
                            b.mat1 = boundaries[j].mat2;
                            b.mat2 = boundaries[j].mat1;
                        }
                        X.push_back(b);
                    }
                }
                // If there are no boundaries in the microlayer...
                if (X.size()==0) {
                    // Then it is a homogeneous layer, and we need to make one...
                    list<bounds>::iterator p,q;
                    for (q=p=Yslice.begin();p!=Yslice.end();++p) {
                        if (p->x>h && q->x<h) {
                            bounds b;
                            b.x=0;
                            b.mat1=p->mat1;
                            b.mat2=p->mat1;
                            X.push_back(b);
                        }
                        q=p;
                    }
                }
                X.sort();
                // Store positions and materials in microlayer
                string imat=X.begin()->mat1;
                for (list<bounds>::iterator q=X.begin();q!=X.end();++q) {
                    if (q->mat1!=imat) 
                        error("Inconsistent material index: " + q->mat1 + "!=" + imat);
                    imat = q->mat2;
                    position[k].push_back(q->x);
                    material[k].push_back(q->mat1);
                }
                if (X.begin()->mat1!=imat) 
                    error("Inconsistent material index: " + X.begin()->mat1 + "!=" + imat);
                position[k].push_back(X.begin()->x+period);
                material[k].push_back(X.begin()->mat1);
                // Increment microlayer number...
                k++;
            }
            prevy=*p;
            // Increment sublayer number...
            ++kk;
        }

        // Grating wants information stored from the top microlayer...
        reverse(position.begin(),position.end());
        reverse(material.begin(),material.end());
        reverse(thickness.begin(),thickness.end());
    }

    COMPLEX 
    Generic_Grating::
    get_complex_value(istream& is) const
    {
        string a = read_pair(is);

        if (is.eof()) error("Premature end of file while reading vertex");

        string b = Evaluator(a,vars).ResultString();

        istringstream c(b);
        COMPLEX d;
        c >> d;
        return d;
    }

    void Generic_Grating::error(const string& message) const
    {
        Model::error("(input file = " + filename + "): " + message);
    }

    COMPLEX Generic_Grating::fourier(int i,int level,int recip)
    {
        SETUP();
        
        Fourier_Component_Calculator z(period,i,recip,Grating::boundary);
        int j;
        for (j=0;j<position[level].size();++j) {
            double x=position[level][j];
            z.enter(x,eps[material[level][j]]);
        }
        COMPLEX result = z.result();
        return result;
    }


    DEFINE_VIRTUAL_MODEL(Grating,Model,"Grating","Generalized description of a grating");
    DEFINE_PARAMETER(Grating,double,period,"Period of the grating [um]","1");
    DEFINE_PARAMETER(Grating,double,boundary,"Graded boundary between materials [um]","0");
    DEFINE_PARAMETER(Grating,dielectric_function,medium_i,"Incident medium","(1,0)");
    DEFINE_PARAMETER(Grating,dielectric_function,medium_t,"Transmitted medium","(4.05,0.05)");

    DEFINE_MODEL(Single_Line_Grating,Grating,"Single_Line_Grating","Model for a single line with adjustable sidewall angles");
    DEFINE_PARAMETER(Single_Line_Grating,dielectric_function,material,"Line material","(4.05,0.05)");
    DEFINE_PARAMETER(Single_Line_Grating,dielectric_function,space,"Space between lines","(1,0)");
    DEFINE_PARAMETER(Single_Line_Grating,double,height,"Height of grating","0.2");
    DEFINE_PARAMETER(Single_Line_Grating,double,topwidth,"Width of top of line [um]","0.2");
    DEFINE_PARAMETER(Single_Line_Grating,double,bottomwidth,"Width of bottom of line [um]","0.2");
    DEFINE_PARAMETER(Single_Line_Grating,double,offset,"Shift of bottom relative to top [um]","0");
    DEFINE_PARAMETER(Single_Line_Grating,int,nlevels,"Number of levels","10");

    DEFINE_MODEL(Corner_Rounded_Grating,Grating,"Corner_Rounded_Grating","A single line grating with corner rounding and a sidewall angle.");
    DEFINE_PARAMETER(Corner_Rounded_Grating,dielectric_function,material,"Line optical properties","(4.05,0.05)");
    DEFINE_PARAMETER(Corner_Rounded_Grating,double,height,"Height of line [um]","0.1");
    DEFINE_PARAMETER(Corner_Rounded_Grating,double,width,"Bottom width of line [um]","0.1");
    DEFINE_PARAMETER(Corner_Rounded_Grating,double,sidewall,"Sidewall angle [deg]","88");
    DEFINE_PARAMETER(Corner_Rounded_Grating,double,radiusb,"Bottom radius [um]","0.010");
    DEFINE_PARAMETER(Corner_Rounded_Grating,double,radiust,"Top radius [um]","0.001");
    DEFINE_PARAMETER(Corner_Rounded_Grating,int,nlevels,"Number of levels","10");

    DEFINE_MODEL(Sinusoidal_Grating,Grating,"Sinusoidal_Grating","A sinusoidal grating");
    DEFINE_PARAMETER(Sinusoidal_Grating,dielectric_function,material,"Line material","(4.05,0.05)");
    DEFINE_PARAMETER(Sinusoidal_Grating,double,amplitude,"Amplitude of sinusoid [um]","0.4");
    DEFINE_PARAMETER(Sinusoidal_Grating,double,base,"Base of sinusoid [um]","0.");
    DEFINE_PARAMETER(Sinusoidal_Grating,int,nlevels,"Number of levels","20");

    DEFINE_MODEL(Triangular_Grating,Grating,"Triangular_Grating","A grating having a triangular profile");
    DEFINE_PARAMETER(Triangular_Grating,dielectric_function,material,"Triangle medium","(4.05,0.05)");
    DEFINE_PARAMETER(Triangular_Grating,double,amplitude,"Amplitude of grating [um]","0.4");
    DEFINE_PARAMETER(Triangular_Grating,double,aspect,"Aspect ratio of triangle (location of apex relative to period)","0.5");
    DEFINE_PARAMETER(Triangular_Grating,int,nlevels,"Number of levels","20");

    DEFINE_MODEL(Generic_Grating,Grating,"Generic_Grating","A grating described by a file specification");
    DEFINE_PARAMETER(Generic_Grating,string,filename,"Filename","");
    DEFINE_PARAMETER(Generic_Grating,string,pstring,"Parameter string","");
    DEFINE_PARAMETER(Generic_Grating,int,nlayers,"Approximate number of levels","20");

    DEFINE_MODEL(Null_Grating,Grating,"Null_Grating","A grating with zero layers.");

}