polefigure.cpp

Go to the documentation of this file.

#include <iostream>
#include <cstdlib>
#include <fstream>
#include <memory>
#include <map>

#include "atomprobe/atomprobe.h"

#include <gsl/gsl_linalg.h>

using namespace std;
using namespace AtomProbe;

void gsl_print_matrix(gsl_matrix *m)
{
    cerr << "[" << endl;
    for (size_t i = 0; i < m->size1; i++)
    {
        for (size_t j = 0; j < m->size2; j++)
        {
            cerr << gsl_matrix_get(m,i,j) << " ";
        }
        cerr << endl;

    }
    cerr << "]" << endl;

}


#include "atomprobe/helper/constants.h"

int main(int argc, char *argv[])
{
    if(! (argc==6  || argc ==7 || argc == 8))
    {
        cerr << "Simulates the projection coordinates in an APT pole figure (simple cubic lattice), using the method of Perry & Brandon [1]. (Make spherical shell and project)" << endl; 
        cerr << endl;
        cerr << "Usage: LATTICENAME SCALE_TO_LATTICE THICKNESS_TO_LATTICE FOCUS FLIGHT_PATH [Forwards_vector] [right_vector|rotation_angle]" << endl;
        cerr << "LATTICENAME : Any of \"Simple\" (simple cubic), FCC or BCC." << endl;
        cerr << "SCALE_TO_LATTICE : This ratio is the density of atoms to project, i.e the number of atoms in the simulation in one direction." << endl;
        cerr << "THICKNESS_TO_LATTICE: Number of lattice slices to take from shell. Moore suggests 0.1" << endl;
        cerr << "FOCUS : The centre projection point for the modified stereographic projection. Sometimes called \"Image compression factor\". Must be >=1" << endl;
        cerr << "FLIGHT_PATH : The units here are arbitrary, and just scale the projection" << endl;

        cerr << "These values are dimensionless, as the projection problem can be expressed as a dimensionlessly (i.e. independently of lattice parameter)" << endl;
        cerr << endl;

        cerr << endl;
        cerr << endl;
        cerr << "Forwards vector : The vector to place at the centre of the pole figure, eg 0,1,0" << endl;
        cerr << "Right vector : The vector to place as the right-most point of the pole figure, eg 1,0,1. If provided, this must be orthogonal to the forwards vector" << endl;
        cerr << " Instead of the right vector, you can provide a rotation angle (degrees). This will rotate the crystal around 0,0,1 (world) after the initial transformation" << endl;
            
        cerr << endl;
        cerr << endl;

        cerr << "[1] A. J. Perry & D. G. Brandon, The bond structure of computer-simulated field-ion images, Philosophical Magazine" << endl;

        cerr << endl;
        cerr << "Example:" << argv[0] << " 10  4 1" << endl;

        cerr << endl;
        cerr << "You can rotate the pole centre by specifying the central atom vector, and the the right hand vector" << endl;
        cerr << "Example:" << argv[0] << " 10  4 (1,1,1) (0,0,1)" << endl;


        cerr << endl;
        cerr << "Output:" << endl;
        cerr << "The output is a 2-column vector with the detector positions. These are in the units given by the FLIGHT_PATH" << endl;
        cerr << "To work out if these would land on any given detector, you must perform this computation externally to this program. All possible visible points from the projection will be shown" << endl;

    }

    //Parse numerical arguments
    //---
    float scaleToLattice,thicknessToLattice,focus,flightPath;
    if(stream_cast(scaleToLattice,argv[2]))
        return 1;
    if(stream_cast(thicknessToLattice,argv[3]))
        return 1;
    if(stream_cast(focus,argv[4]))
        return 1;
    if(stream_cast(flightPath,argv[5]))
        return 1;
    //---


    //Check if user specified orientation vectors
    //----
    bool specifiedForwards=false;
    if(argc == 7)
        specifiedForwards=true;
    bool specifiedRight=false;
    if(argc == 8)
        specifiedRight=true;
    //----

    //Check argument validity
    //--
    if(scaleToLattice <= 1e-5)
    {
        cerr << "Scale/Lattice ratio too small!" << endl;
        return 1;
    }

    if(thicknessToLattice <=1e-5)
    {
        cerr << "Thickness/lattice ratio too small" << endl;
        return 1;
    }


    if(focus <=-1)
    {
        cerr << "Projection focus must be greater than -1" << endl;
        return 1;
    }

    if(flightPath < 0)
    {
        cerr << "Flight path must be positive" << endl;
        return 1;
    }
    //--


    //See if we want to post-rotate the crystal
    float postRotateAngle=0.0f;
    //Parse orientation vectors, if supplied
    //----
    bool rotate=false;
    Point3D forwardsVec,rightVec;
    if(specifiedForwards || specifiedRight)
    {
        if(!forwardsVec.parse(argv[6]))
        {
            cerr << "Error parsing  forwards vector (centre of figure)" << endl;
            return 1;
        }
        forwardsVec.normalise();
    
        if(specifiedRight)
        {
            if(!rightVec.parse(argv[7]))
            {
                if(stream_cast(postRotateAngle,argv[7]))
                {
                    cerr << "Error parsing  right vector/rotation angle " << endl;
                    return 1;
                }
                else
                {
                    //Convert from deg to rad
                    postRotateAngle*=M_PI/180.0f;
                    rightVec=Point3D(0,0,1).crossProd(forwardsVec);
                    cerr << "Right vector is :" << rightVec << endl;
                }
            }
        }
        else
        {
            //Compute the right vector as normal to 0,0,1 and our new forwards vector
            rightVec=Point3D(0,0,1).crossProd(forwardsVec);
            cerr << "Right vector is :" << rightVec << endl;
        }

        rightVec.normalise();
        
        if(fabs(forwardsVec.sqrMag() -1.0f) > 1e-3 ||
            fabs(rightVec.sqrMag() -1.0f) > 1e-3)
        {
            cerr << "Input vectors could not be normalised." << endl;
            return 1;
        }

        if(fabs(forwardsVec.dotProd(rightVec)) > 1e-3)
        {
            cerr << "Requested pole figure vectors are not orthogonal. " << endl;
            return 1;
        }

        rotate=true;
    }
    else
    {
        forwardsVec=Point3D(0,0,1);
        rightVec=Point3D(1,0,0);
    }
    //----


    //Convert name of lattice to type identifier
    //--
    enum
    {
        LATTICE_SIMPLE_CUBIC,
        LATTICE_FCC,
        LATTICE_BCC,
    };

    map<string,unsigned int> latticeNameMap;
    latticeNameMap["simple"] = LATTICE_SIMPLE_CUBIC;
    latticeNameMap["fcc"] = LATTICE_FCC;
    latticeNameMap["bcc"] = LATTICE_BCC;


    string latticeName=argv[1];
    latticeName=lowercase(latticeName);
    
    unsigned int latticeType;
    auto it =latticeNameMap.find(latticeName); 
    if(it== latticeNameMap.end())
    {
        cerr << "Unsupported or unknown lattice:" << latticeName << endl;
        return 1;
    }

    latticeType=it->second;


    
    //Generate a lattice in the upper quadrant  
    Point3D p(scaleToLattice,scaleToLattice,scaleToLattice);
    
    //Select the crystal lattice to generate
    CrystalGen *gen;
    switch(latticeType)
    {
        case LATTICE_SIMPLE_CUBIC:
            gen=(new SimpleCubicGen(1,1,p));
            break;
        case LATTICE_FCC:
        {
            const float massToCharge[3]={1,1,1};
            gen=(new FaceCentredCubicGen(1,1,massToCharge,p));
            break;
        }
        case LATTICE_BCC:
        {
            gen=(new BodyCentredCubicGen(1,1,2,p));
            break;
        }
        default:
        {
            cerr << "Unknown lattice. Aborting" << endl;
            return 1;
        }
    }

    gen->generateLattice();
    //Copy the lattice into all 8 quadrants
    gen->mirrorOut();

    vector<Point3D> h;
    gen->extractPositions(h);
    cerr << "Full lattice contains :" << h.size() << "pts" << endl;
    
    delete gen;

    //Perform rotation, if required
    //==
    if(rotate)
    {

        //Compute the rotation matrix between the two vectors
        gsl_matrix *rMat = gsl_matrix_alloc(3,3);
        computeRotationMatrix(Point3D(0,0,1),Point3D(1,0,0),forwardsVec,rightVec,rMat);

    
        //Compute the inverse matrix to take us from {forwardsvec,rightvec} -> {(0,0,1),(1,0,0)}
        // The inverse of an orthonomal matrix is transpose
        gsl_matrix_transpose(rMat);

        for(unsigned int ui=0;ui<h.size();ui++)
            h[ui].transform3x3(rMat);

        //User requested post-rotation  to bring the pole figure into alignment
        if(postRotateAngle!=0.0f)
        {
            float sR = sin(postRotateAngle);
            float cR = cos(postRotateAngle);
#pragma omp parallel for
            for(unsigned int ui=0;ui<h.size();ui++)
            {
                //Perform a 2D rotation of H
                Point3D pt;
                pt = h[ui];
                pt[0] = h[ui][0]*cR - h[ui][1]*sR;
                pt[1] = h[ui][0]*sR + h[ui][1]*cR;
                h[ui] = pt;
            }

            //compute the above rotation as a matrix
            gsl_matrix *m=gsl_matrix_alloc(3,3);
            gsl_matrix_set_identity(m);
            gsl_matrix_set(m,0,0,cR);
            gsl_matrix_set(m,0,1,-sR);
            gsl_matrix_set(m,1,0,sR);
            gsl_matrix_set(m,0,1,cR);

            //transpose to invert the operaiton
            gsl_matrix_transpose(m);
            
            //multiply this by rMat 
            gsl_matrix *tmp=gsl_matrix_alloc(3,3);
            gsl_matrix_mult(rMat,m,tmp);

            std::swap(tmp,rMat);
            gsl_matrix_free(tmp);
            
        }

        //Transpose to obtain reverse transformation matrix
        gsl_matrix_transpose(rMat);

        //Print the reverse matrix
        cerr << "Reverse transformation (World -> crystal)" << endl;
        cerr << "-----------------" << endl;
        gsl_print_matrix(rMat);
        cerr << "-----------------" << endl;
        
        gsl_matrix_free(rMat);
    }

    //==
    //Keep only the upper plane
    {

    vector<Point3D> quadPts; 
    for(unsigned int ui=0;ui<h.size();ui++)
    {
        if(h[ui][2] >=0)
            quadPts.push_back(h[ui]);
    }

    h.swap(quadPts);
    }


    //Compute which points are to stay within the slice.
    // discard the rest 
    // FIXME: This would be better written as lambda function passed to the generator
    //      to reduce memory requirements
    vector<Point3D> pts;
    float deltaMin=(scaleToLattice-thicknessToLattice)*(scaleToLattice-thicknessToLattice);
    float deltaMax=(scaleToLattice)*(scaleToLattice);
    for(size_t ui=0;ui<h.size();ui++)
    {
        float rSqr;
        rSqr = ( h[ui].sqrMag());
        
        if( rSqr > deltaMin && rSqr < deltaMax)
            pts.push_back(h[ui]);
             
    }
    cerr << "Sliced out :" << pts.size() << " pts" << endl; 

    if(!pts.size())
    {
        cerr << "No points outputted. Consider increasing your lattice density (scale/lattice) or slice thickness" << endl;
        return 1;
    }


    savePosFile(pts,1,"sliced-lattice.pos");
    
    ModifiedFocusSphericProjection mSp(focus);  
    vector<pair<float,float> > pXY;
    //The FOV reported by maxFOV is the full angle, not the half-angle
    float maxFOV = mSp.getMaxFOV()/2.0f;
    for(size_t ui=0;ui<pts.size() ;ui++)
    {
        //Convert from the spherical slice to spherical coordinates
        float eta,phi,r;
        pts[ui].getISOSpherical(r,eta,phi);

        //Convert from spherical coordinates to projection coords
        float theta;
        theta=mSp.etaToTheta(eta);
        //project these coordinates using the above projection
        float pX,pY;
        if( eta < maxFOV && mSp.toPlanar(theta,phi,pX,pY))
        {
            pXY.push_back(make_pair(pX,pY));
        }
    }

    //Conversion factor to scale up to detector coordinates 
    // this arises from similar triangles, where one triangle
    // is the m+1 (x-axis) and projected coord (p)
    float flightFactor=flightPath/(focus+1);
    

    ofstream fProj("projection.txt");
    if(!fProj)
    {
        cerr << "Failed to open output file, aborting" << endl;
        return 1;
    }


    //Print arguments used to call program to file
    fProj << "# Arguments used: ";
    for(unsigned int ui=1;ui<argc;ui++)
        fProj << argv[ui] << " ";
    fProj << endl;

    //record the projected XY coords, scaling up by the flight factor
    for(size_t ui=0;ui<pXY.size();ui++)
    {
        fProj << pXY[ui].first*flightFactor << " , " << pXY[ui].second*flightFactor << endl;
    }

    cerr << "Wrote detector coordinates to projection.txt" << endl;

    cerr << "Max FOV:" << mSp.getMaxFOV()*180.0/M_PI << endl;

    



    return 0;
}