/*******************************************************************************
*
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: MM_c
*
* %Identification
* Written by: Erik B Knudsen and Desiree D. M. Ferreira
* Date: Feb. 2016, Feb. 2017
* Modified by: Søren Jeppesen 
* Version: 1.1
* Release: McXtrace 1.2
* Origin: DTU Physics, DTU Space
*
* Single Pore as part of the Silicon Pore Optics (SPO) as envisioned for the ATHENA+ space telescope.
*
* %Description
* A single pore is simulated, which may have thick walls. The top and bottom are curved cylindrically
* azimuthally, and according to the Wolter I optic lengthwise (sagitally). A parameter specifies
* whether this is hyperbolic or parabolic.
* The azimuthal curvature is defined by the parameter radius. This refers to the center of the pore. I.e the top
* and bottom plates have radius of curvature <radius+yheight/2> and <radius-yheight/2> respectively.
*
* To intersect the Wolter I plates we take advatage of the azimuthal symmetry and only consider the radial component
* of the photon's wavevector.
*
* %Parameters
* Input parameters:
* radius_m: [m] Ring radius of the upper (reflecting) plate of the pore at the optic centre.
* yheight:  [m] Height of the pore.
* xwidth:   [m] Width of the pore.
* chamferwidth:  [m] Width of side walls.
* gap:      [m] Gap between the plate and the intersection plane with the hyperbolic section. (currently ignored)
* Z0:       [m] Distance between optics centre plane and focal spot (essentially focal length).
* mirror_reflec: [ ] Data file containing reflectivities of the reflector surface (TOP).
* side_reflec:   [ ]   Data file containing reflectivities of the side walls (LEFT and RIGHT).
* bottom_reflec: [ ]  Data file containing reflectivities of the bottom surface (BOTTOM).
* R_d:      [ ] Default reflectivity value to use if no reflectivity file is given. Useful f.i. is one surface is reflecting and the others absorbing.
* primary:  [ ] If non-zero, the pore is considered a primary reflector, and extends towards negative z. I.e. the entry plane is behind the z=0-plane. If zero, the pore is considered secondary and extends from the z=0-plane and towards positive z.
* dalpha:   [deg] Offset to the alpha angle computed from the focal length. Useful for targeting the modified conical geometry (currently ignored).
* waviness: [rad] Waviness of the reflecting surface. The slope error is assumed to be uniformly distributed in the interval [-waviness:waviness].
* longw:    [ ] If non-zero, waviness is 1D and along the pore axis.
* %End
*******************************************************************************/

DEFINE COMPONENT MM_c
SETTING PARAMETERS (pore_th, pore_width, int ring_nr,radius_m, Z0, xwidth, pore_height, gap=0, chamfer_width=0, length=0, string mirror_reflec="", string bottom_reflec="", string side_reflec="", string size_file="", string non_specular_file="", R_d=1, primary=1, dalpha=0, waviness=0, longw=0)

SHARE
%{
#ifndef MCSPO_INTERSECT_CONE
#define MCSPO_INTERSECT_CONE 1

    int intersect_cone(double *l0, double x, double y, double z, double kx, double ky, double kz, double alpha, double radius, double *nx, double *ny, double *nz){
        double kxn=kx,kyn=ky,kzn=kz;
        NORM(kxn,kyn,kzn);
        double c=tan(alpha);
        double z0=radius/c;
        double c2=c*c;
        double A,B,C;
        A=kxn*kxn + kyn*kyn - c2*kzn*kzn;
        B=2*(kxn*x + kyn*y  - c2*kzn*(z-z0));
        C=x*x + y*y - c2*(z-z0)*(z-z0);

        int status;
        double l1;
        if ( (status=solve_2nd_order(l0,&l1,A,B,C))==0 ){
            /*note that if l1->NULL only the smallest positive solution is returned*/
            //fprintf(stderr,"Error(%s): No solution to second order eq.\n","MM_c");
            return status;
        }
        /*compute normal vector here*/
        x+=kxn* (*l0);
        y+=kyn* (*l0);
        z+=kzn* (*l0);

        double vn=sqrt(x*x+y*y);
        *nx=x/vn;
        *ny=y/vn;

        *nz=1;

        *nx *= cos(alpha);
        *ny *= cos(alpha);
        *nz *= sin(alpha);

        return status;
    }
#endif

#ifndef PROP_Z
#define PROP_Z(zz)\
    mcPROP_Z(zz)
#endif

#ifndef mcPROP_Z
#define mcPROP_Z(zz) \
  do { \
    MCNUM mc_dl,mc_k; \
    if(kz == 0) { ABSORB; }; \
    mc_k=sqrt(scalar_prod(kx,ky,kz,kx,ky,kz)); \
    mc_dl= ((zz)-z) * mc_k / kz; \
    if(mc_dl<0 && mcallowbackprop==0) { ABSORB; };\
    PROP_DL(mc_dl); \
  } while(0)
#endif

#ifndef DISTRIBUTIONS
#define DISTRIBUTIONS
    double laplaceDistribution(double mean, double width, _class_particle* _particle){
        double sample = rand01();
        if(sample <= 0.5){
            return mean + width*log(2*sample);
        } else {
            return mean - width*log(2-2*sample);
        }
    }

    double lorentzDistribution(double mean, double width, _class_particle* _particle){
        double sample = rand01();
        return mean + width*tan(M_PI*(sample-0.5));
    }
#endif
%}

DECLARE
%{
    double nExit[3];
    double wExit[3];
    double nEntry[3];
    double wEntry[3];
    double nTop[3];
    double nBottom[3];
    double nRight[3];
    double wRight[3];
    double nLeft[3];
    double wLeft[3];
    double e_min[3];
    double e_step[3];
    double e_max[3];
    double theta_min[3];
    double theta_step[3];
    double theta_max[3];

    double e_min_nonspec;
    double e_max_nonspec;
    double e_step_nonspec;
    double theta_min_nonspec;
    double theta_max_nonspec;
    double theta_step_nonspec;


    double zentry;
    double zexit;
    double *zentry_vec;
    double *zexit_vec;
    double *radius_1_vec;
    double *radius_2_vec;

    t_Table reflec_top_table;
    t_Table reflec_bottom_table;
    t_Table reflec_side_table;
    t_Table size_table;
    // array which stores the pore walls that the ray reflects onto
    unsigned char reflections[16];
    
    // index of the latest reflection to occur
    char ref_index;

    t_Table non_specular_table;
%}

INITIALIZE
%{
    #define MAX_TRACKED_REFLECTIONS 16
    char *filenames[]={mirror_reflec,bottom_reflec,side_reflec};
    t_Table *ref_tables[]={&reflec_top_table,&reflec_bottom_table,&reflec_side_table};
    int i;

    /*read data from files into tables using read_table-lib*/
    for (i=0;i<3;i++){
        char *reflec=filenames[i];
        t_Table *tp=ref_tables[i];
        if (reflec && strlen(reflec)) {
            char **header_parsed;

            /* read 1st block data from file into tp */
            if (Table_Read(tp, reflec, 1) <= 0)
            {
                exit(fprintf(stderr,"Error(%s): cannot read file %s\n",NAME_CURRENT_COMP, reflec));
            }
            header_parsed = Table_ParseHeader(tp->header,
                    "e_min=","e_max=","e_step=","theta_min=","theta_max=","theta_step=",NULL);
            if (header_parsed[0] && header_parsed[1] && header_parsed[2] &&
                    header_parsed[3] && header_parsed[4] && header_parsed[5])
            {
                e_min[i]=strtod(header_parsed[0],NULL);
                e_max[i]=strtod(header_parsed[1],NULL);
                e_step[i]=strtod(header_parsed[2],NULL);
                theta_min[i]=strtod(header_parsed[3],NULL);
                theta_max[i]=strtod(header_parsed[4],NULL);
                theta_step[i]=strtod(header_parsed[5],NULL);
            } else {
                exit(fprintf(stderr,"Error (%s): wrong/missing header line(s) in file %s\n", NAME_CURRENT_COMP, reflec));
            }
            if (!((int)(e_max[i]-e_min[i]) == (int)((tp->rows-1)*e_step[i])))
            {
                exit(fprintf(stderr,"Error (%s): e_step does not match e_min and e_max in file %s\n",NAME_CURRENT_COMP, reflec));
            }
            if (!((int)(theta_max[i]-theta_min[i]) == (int)((tp->columns-1)*theta_step[i])))
            {
                exit(fprintf(stderr,"Error (%s): theta_step does not match theta_min and theta_max in file %s\n",NAME_CURRENT_COMP, reflec));
            }
        }else{
            /*mark the table as unread by setting "rows" to -1
              This will trigger the default reflectivity.*/
            tp->rows=-1;
        }
    }

    /* Read table with mirror module parameters*/

    if(size_file){
    	if (Table_Read(&(size_table), size_file, ring_nr) <=0)
    	  exit(fprintf(stderr, "Error (%s): Could not read %s. Aborting.\n",NAME_CURRENT_COMP, size_file));
    }

    zentry_vec=calloc(size_table.rows,sizeof(double));
    zexit_vec=calloc(size_table.rows,sizeof(double));
    radius_1_vec=calloc(size_table.rows,sizeof(double));
    radius_2_vec=calloc(size_table.rows,sizeof(double));
    /*all the plates in 1 ring are identical - i.e. we can set the plate xwidth once and for all.*/
    xwidth=Table_Index(size_table,0,6);

    double alpha,thetap,thetah,P,d,e,C0,Z;
    double radius_m, radius_1, radius_2;
    /*There are in general 68 reflecting planes*/
    for (i=0;i<size_table.rows;i++){

        radius_m = Table_Index(size_table,i,4);

        /* compute some pore parameters*/
        alpha=0.25*atan(radius_m/Z0);

        double D,Z0p;
        if (primary){
	    Z0p=radius_m/tan(alpha);
            length=Table_Index(size_table,i,1);
            D=sqrt(radius_m*radius_m + Z0p*Z0p);
            zentry=Z0p*(1-(D+length)/D);
            zexit=0;
            radius_1=(D+length)/D*radius_m;
            radius_2=radius_m;
        }else{
            alpha*=3.0;
            Z0p=radius_m/tan(alpha);
            length=Table_Index(size_table,i,2);
            D=sqrt(radius_m*radius_m + Z0p*Z0p);
            zentry=0;
            zexit=Z0p*(1-(D-length)/D);
            radius_1=radius_m;
            radius_2=(D-length)/D*radius_m;
        }
        zentry_vec[i] = zentry;
        zexit_vec[i] = zexit;
        radius_1_vec[i] = radius_1;
        radius_2_vec[i] = radius_2;
    }

    /*find minimum zentry - i.e. the "first" entry plane.*/
    zentry = zentry_vec[0];
    for(i=0;i<size_table.rows-1;i++){
        if(zentry_vec[i] < zentry){
            zentry = zentry_vec[i];
        }
    }

    zexit = zexit_vec[0];
    for(i=0;i<size_table.rows-1;i++){
        if(zexit_vec[i] > zexit){
            zexit = zexit_vec[i];
        }
    }

    nEntry[0]=0;
    nEntry[1]=0;
    nEntry[2]=-1;
    wEntry[0]=wEntry[1]=0;wEntry[2]=zentry;

    nExit[0]=0;
    nExit[1]=0;
    nExit[2]=1;
    wExit[0]=wExit[1]=0;wExit[2]=zexit;
    if(non_specular_file && strlen(non_specular_file)){
        if (Table_Read(&(non_specular_table), non_specular_file, 0) <= 0)
          exit(fprintf(stderr, "Error (%s): Could not read %s. Aborting.\n",NAME_CURRENT_COMP, non_specular_file));
        char** header_parsed = Table_ParseHeader(non_specular_table.header,
            "e_min=","e_max=","e_step=","theta_min=","theta_max=","theta_step=",NULL);
            if (header_parsed[0] && header_parsed[1] && header_parsed[2] &&
                header_parsed[3] && header_parsed[4] && header_parsed[5])
                {
                    e_min_nonspec=strtod(header_parsed[0],NULL);
                    e_max_nonspec=strtod(header_parsed[1],NULL);
                    e_step_nonspec=strtod(header_parsed[2],NULL);
                    theta_min_nonspec=strtod(header_parsed[3],NULL);
                    theta_max_nonspec=strtod(header_parsed[4],NULL);
                    theta_step_nonspec=strtod(header_parsed[5],NULL);
                } else {
                    exit(fprintf(stderr,"Error (%s): wrong/missing header line(s) in file %s\n", NAME_CURRENT_COMP, non_specular_file));
                }
                if (!((int)(e_max_nonspec-e_min_nonspec) == (int)((non_specular_table.rows-1)*e_step_nonspec)))
                {
                    exit(fprintf(stderr,"Error (%s): e_step does not match e_min and e_max in file %s\n",NAME_CURRENT_COMP, non_specular_file));
                }
                if (!((int)(theta_max_nonspec-theta_min_nonspec) == (int)((non_specular_table.columns-1)*theta_step_nonspec)))
                {
                    exit(fprintf(stderr,"Error (%s): theta_step does not match theta_min and theta_max in file %s\n",NAME_CURRENT_COMP, non_specular_file));
                }
    }
%}

TRACE
%{
    ref_index = -1;
    double r_entry_min,r_entry_max,r2;
    int hit=0;
    int hit_chamfer=0;
    double psi_max,psi_min,psi;
    double radius_m;

    /*assuming the table to be sorted*/
    ALLOW_BACKPROP;
    PROP_Z(zentry);
    r_entry_min = radius_1_vec[0];
    r_entry_max = radius_1_vec[size_table.rows-1];

    r2 = (x*x + y*y);

    //bounds of the module as a whole
    psi_min=-xwidth*0.5/(r_entry_min-pore_height);
    psi_max= xwidth*0.5/(r_entry_min-pore_height);
    psi=atan2(x,y);

    hit = ( ( r2 > ( r_entry_min-pore_height )*( r_entry_min-pore_height ) ) && ( r2 < ( r_entry_max + pore_th )*( r_entry_max + pore_th ) )
            && (psi>psi_min && psi<psi_max) );
    if (hit){
        /*we are within range for entering the the mirror module - we might still miss due to beam divergence though*/
        hit=0;
        int ii=0;
        int hit_chamfer=0;
        int upside_down_plate;
        double local_pore_height;
        double radius_entry;

        while( !hit && ii<size_table.rows){
            ALLOW_BACKPROP;
            PROP_Z(zentry_vec[ii]);

            radius_m = Table_Index(size_table,ii,4);
            radius_entry = radius_1_vec[ii];

            upside_down_plate = Table_Index(size_table, ii, 8);
            if(upside_down_plate){
                local_pore_height = 2*pore_height;
            } else {
                local_pore_height = pore_height;
            }

            hit= ( ( x*x + y*y < radius_entry*radius_entry ) && ( x*x + y*y >(radius_entry-local_pore_height)*(radius_entry-local_pore_height) ) ) ;
            ii++;
        }
        /*if we have missed all plates, terminate the ray*/
        if (!hit){
            ABSORB;
        }

        double alpha=0.25*atan(radius_m/Z0);
        if(!primary){
            alpha *= 3;
        }

        /*figure out which pore we hit*/
        int jj=0;
        // Width of the plate in pores
        double width_pores = (xwidth-chamfer_width)/(pore_width+chamfer_width);

        // photon psi in pores
        double psi_pores = (psi - (-xwidth*0.5)/(radius_entry-local_pore_height)) *(radius_entry-local_pore_height)/ ( radius_entry+chamfer_width);
        char hit_final_pore = 0;
        if(psi_pores >= 0 && psi_pores < width_pores){
            hit = 1 & hit;
            jj = (int) floor(psi_pores);
            //If the pore that is hit is the last pore on the plate, modify the
            //pore width and chamfer width so it corresponds to the final, larger
            //pore.
            if(jj >= floor(width_pores) - 1){
                jj = floor(width_pores) - 1;
                hit_final_pore = 1;
            }
        } else {
            hit = 0;
        }
        /*check for side_wall (chamfer) hit (obscuration)*/
        hit_chamfer=0;

        if(hit_final_pore && psi > xwidth/2/(radius_entry-local_pore_height) - chamfer_width){
            hit_chamfer=1;
            hit=0;
            ABSORB;
        } else if ( (psi - (-xwidth/2.0 + jj*(pore_width+chamfer_width))/(radius_entry-local_pore_height)) <chamfer_width/(radius_entry-local_pore_height) ){
            hit_chamfer=1;
            hit=0;
            ABSORB;
        }

        int coatedPlate = Table_Index(size_table, ii, 7);

        enum {LEFT, RIGHT, TOP, BOTTOM, EXIT, NONE} wall;
        t_Table *reflec_table=NULL;
        double R;

        if(hit){
            SCATTER;
            int exit=0;
            int intersections[5];
            int i_small;
            double l[5];
            double l_small;

            double nx,ny,nz;

            psi_min=( -xwidth/2.0 + jj*(pore_width+chamfer_width) + chamfer_width)/(radius_entry-local_pore_height);
            if (hit_final_pore){
                psi_max = (xwidth/2.0-chamfer_width)/(radius_entry-local_pore_height);
            } else {
                psi_max=( -xwidth/2.0 + (jj+1)*(pore_width+chamfer_width))/(radius_entry-local_pore_height);
            }

            /*TODO side wall planes*/
            nLeft[0]=cos(psi_min); nLeft[1]=-sin(psi_min); nLeft[2]=0;
            wLeft[0]=radius_entry*sin(psi_min); wLeft[1]=radius_entry*cos(psi_min); wLeft[2]=0;

            nRight[0]=cos(psi_max); nRight[1]=-sin(psi_max); nRight[2]=0;
            wRight[0]=radius_entry*sin(psi_max); wRight[1]=radius_entry*cos(psi_max); wRight[2]=0;

            while (!exit){
                l_small=DBL_MAX;
                wall=NONE;
                double nx,ny,nz;
                double wx,wy,wz;
                int prm_idx;/*index indicating which table parameter set to choose*/

                /*left wall*/
                intersections[LEFT]=plane_intersect(l+LEFT,x,y,z,kx,ky,kz,nLeft[0],nLeft[1],nLeft[2],wLeft[0],wLeft[1],wLeft[2]);
                if (intersections[LEFT] && l[LEFT]>DBL_EPSILON && l[LEFT]<l_small) {l_small=l[LEFT];i_small=intersections[LEFT];wall=LEFT;}
                /*right wall*/
                intersections[RIGHT]=plane_intersect(l+RIGHT,x,y,z,kx,ky,kz,nRight[0],nRight[1],nRight[2],wRight[0],wRight[1],wRight[2]);
                if (intersections[RIGHT] && l[RIGHT]>DBL_EPSILON && l[RIGHT]<l_small) {l_small=l[RIGHT];i_small=intersections[RIGHT];wall=RIGHT;}
                /*exit plane*/
                intersections[EXIT]=plane_intersect(l+EXIT,x,y,z,kx,ky,kz,nExit[0],nExit[1],nExit[2],wExit[0],wExit[1],wExit[2]);
                if (intersections[EXIT] && l[EXIT]>DBL_EPSILON && l[EXIT]<l_small) {l_small=l[EXIT];i_small=intersections[EXIT];wall=EXIT;}
                /*top surface - the real reflecting surface*/
                intersections[TOP]=intersect_cone((l+TOP),x,y,z,kx,ky,kz,alpha,radius_m,&(nTop[0]),&(nTop[1]),&(nTop[2]));
                if (intersections[TOP] && l[TOP]>DBL_EPSILON && l[TOP]<l_small) {l_small=l[TOP];i_small=intersections[TOP];wall=TOP;}
                /*bottom surface*/
                intersections[BOTTOM]=intersect_cone((l+BOTTOM),x,y,z,kx,ky,kz,alpha,radius_m-local_pore_height,&(nBottom[0]),&(nBottom[1]),&(nBottom[2]));
                if (intersections[BOTTOM] && l[BOTTOM]>DBL_EPSILON && l[BOTTOM]<l_small) {l_small=l[BOTTOM];i_small=intersections[BOTTOM];wall=BOTTOM;}

                /*sort intersections to find the smallest positive one*/
                switch (wall){
                    case LEFT:
                        /*handle left wall "reflection"*/
                        reflec_table=&reflec_side_table;
                        nx=nLeft[0];ny=nLeft[1];nz=nLeft[2];
                        prm_idx=2;
                        break;
                    case RIGHT:
                        /*handle right wall "reflection"*/
                        reflec_table=&reflec_side_table;
                        nx=nRight[0];ny=nRight[1];nz=nRight[2];
                        prm_idx=2;
                        break;
                    case TOP:
                        /*handle top wall reflection*/
                        if(coatedPlate){
                            reflec_table=&reflec_top_table;
                            prm_idx=0;
                        } else {
                            reflec_table=&reflec_bottom_table;
                            prm_idx=1;
                        }
                        nx=nTop[0];ny=nTop[1];nz=nTop[2];
                        break;
                    case BOTTOM:
                        /*handle bottom wall "reflection"*/
                        reflec_table=&reflec_bottom_table;
                        nx=nBottom[0];ny=nBottom[1];nz=nBottom[2];
                        prm_idx=1;
                        break;
                    case EXIT:
                        /*photon will exit pore*/
                        exit=1;
                        break;
                }

                ref_index++;

                if(exit){
                    break;
                } else if(ref_index < MAX_TRACKED_REFLECTIONS){
                  reflections[ref_index] = wall;
                }
                PROP_DL(l_small);

                double kix=kx,kiy=ky,kiz=kz;
                double k=sqrt(kx*kx+ ky*ky + kz*kz);
                double e=K2E*k;
                double s=scalar_prod(kx,ky,kz,nx,ny,nz);
                double theta=RAD2DEG*(M_PI_2-acos(s/k)); /*pi_2 since theta is supposed to be the grazing angle*/

                /*if we have waviness alter the normal vector slightly*/
                if(waviness!=0){
                    /*assuming theta to be small we might disregard atan*/
                    if(longw){
                        double dtheta;
                        if(&non_specular_table!=NULL){
                            double lw;
                            lw = Table_Value2d(non_specular_table, (e-e_min[2])/e_step[2], (theta-theta_min[2])/theta_step[2]);
                            dtheta=lorentzDistribution(0, lw, _particle);
                            printf("lorentz dtheta %e\n");
                        } else if(theta<waviness){
                            dtheta=rand01()*(theta+waviness)-theta;
                        } else {
                            dtheta=randpm1()*waviness;
                        }
                        double tx,ty,tz;
                        vec_prod(tx,ty,tz,0,0,1,nx,ny,nz);
                        rotate(nx,ny,nz, nx,ny,nz, dtheta, tx,ty,tz);
                    }else{
                        /*waviness is also transversal but isotropic*/
                        double radius;
                        if(theta<waviness){
                            radius=atan(waviness);
                            randvec_target_circle(&nx,&ny,&nz,NULL,nx,ny,nx,radius);
                        }else{
                            radius=(atan(theta)+atan(waviness))/2.0;
                            randvec_target_circle(&nx,&ny,&nz,NULL,nx,ny,nx+radius-atan(theta),radius);
                        }
                        NORM(nx,ny,nz);
                    }/*reflect the photon through the surface normal*/
                    /*recompute theta*/
                    theta=RAD2DEG*0.5*acos(scalar_prod(kx,ky,kz,kix,kiy,kiz)/k/k);
                }
                /*reflect the photon through the surface normal*/
                if(s!=0){
                    kx-=2*s*nx;
                    ky-=2*s*ny;
                    kz-=2*s*nz;
                }

                if(&non_specular_table!=NULL){
                    double lw = Table_Value2d(non_specular_table, (e-e_min_nonspec)/e_step_nonspec, (theta-theta_min_nonspec)/theta_step_nonspec);
                    double dtheta=lorentzDistribution(0, lw, _particle)*M_PI/180;
                    double tx,ty,tz;
                    vec_prod(tx,ty,tz,0,0,1,nx,ny,nz);
                    rotate(kx, ky, kz, kx, ky, kz, dtheta, tx,ty,tz);
                }

                SCATTER;

                if(reflec_table==NULL || reflec_table->rows==-1){
                    R=R_d;
                }else{
                    R=Table_Value2d(*reflec_table,(e-e_min[prm_idx])/e_step[prm_idx], (theta-theta_min[prm_idx])/theta_step[prm_idx]);
                }
                p*=R;
            }
        }else if (hit_chamfer){
            ABSORB;
        }else{
            /*no hit*/
            ABSORB;
        }
    }else{
        ABSORB;
    }
%}

FINALLY
%{
    free(zentry_vec);
%}

MCDISPLAY
%{
    magnify("");

    int k,j;
    double dth,theta,t0,t1,inner_m,inner_p;
    double r_p[3],r_m[3],w[3],yh[2];
    const int N=20;

    r_p[0] = Table_Index(size_table,size_table.rows,5);
    r_p[1] = Table_Index(size_table,size_table.rows/2,5);
    r_p[2] = Table_Index(size_table,0,5);
    yh[0]=r_p[0]-r_p[2];
    yh[1]=pore_height;

    r_m[0] = Table_Index(size_table,size_table.rows,4);
    r_m[1] = Table_Index(size_table,size_table.rows/2,4);
    r_m[2] = Table_Index(size_table,0,4);

    w[0]=w[2]=xwidth;w[1]=pore_width;

    for (j=0;j<2;j++){
        theta=w[j]/2.0/r_m[j];
        inner_m=r_m[j]-yh[j];
        inner_p=r_p[j]-yh[j];

        line(0,0+r_m[j],0, 0,r_p[j]-r_m[j]+r_m[j],zentry); /*this extra line indicates the reflecting surface*/
        line( sin(theta)*r_m[j], (cos(theta)-1)*r_m[j]+r_m[j], 0,  sin(theta)*r_p[j], cos(theta)*r_p[j]-r_m[j]+r_m[j], zentry);
        line(-sin(theta)*r_m[j], (cos(theta)-1)*r_m[j]+r_m[j], 0, -sin(theta)*r_p[j], cos(theta)*r_p[j]-r_m[j]+r_m[j], zentry);

        line( sin(theta)*inner_m, cos(theta)*inner_m-r_m[j]+r_m[j], 0,  sin(theta)*inner_p, cos(theta)*inner_p-r_m[j]+r_m[j], zentry);
        line(-sin(theta)*inner_m, cos(theta)*inner_m-r_m[j]+r_m[j], 0, -sin(theta)*inner_p, cos(theta)*inner_p-r_m[j]+r_m[j], zentry);

        line( sin(theta)*r_m[j], (cos(theta)-1)*r_m[j]+r_m[j], 0,  sin(theta)*inner_m, cos(theta)*inner_m-r_m[j]+r_m[j], 0);
        line(-sin(theta)*r_m[j], (cos(theta)-1)*r_m[j]+r_m[j], 0, -sin(theta)*inner_m, cos(theta)*inner_m-r_m[j]+r_m[j], 0);
        line( sin(theta)*r_p[j], cos(theta)*r_p[j]-r_m[j]+r_m[j], zentry,  sin(theta)*inner_p, cos(theta)*inner_p-r_m[j]+r_m[j], zentry);
        line(-sin(theta)*r_p[j], cos(theta)*r_p[j]-r_m[j]+r_m[j], zentry, -sin(theta)*inner_p, cos(theta)*inner_p-r_m[j]+r_m[j], zentry);

        dth=2*theta/N;
        for (k=1;k<N+1;k++){
            t0=-theta+(k-1)*dth;
            t1=-theta+k*dth;
            line( sin(t0)*r_m[j], cos(t0)*r_m[j]-r_m[j]+r_m[j], 0, sin(t1)*r_m[j], cos(t1)*r_m[j]-r_m[j]+r_m[j], 0);
            line( sin(t0)*inner_m, cos(t0)*inner_m-r_m[j]+r_m[j], 0, sin(t1)*inner_m, cos(t1)*inner_m-r_m[j]+r_m[j], 0);

            line( sin(t0)*r_p[j], cos(t0)*r_p[j]-r_m[j]+r_m[j], zentry, sin(t1)*r_p[j], cos(t1)*r_p[j]-r_m[j]+r_m[j], zentry);
            line( sin(t0)*inner_p, cos(t0)*inner_p-r_m[j]+r_m[j], zentry, sin(t1)*inner_p, cos(t1)*inner_p-r_m[j]+r_m[j], zentry);
        }
    }
%}

END