temCalcsDemo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math as ma\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle\n",
    "import re\n",
    "from fractions import Fraction\n",
    "import math\n",
    "\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import temCalcs as tc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "steel = tc.addStructure(name=\"steel\", a=3.66, b=3.66, c=3.66, alpha=90, beta=90, gamma=90)\n",
    "dummy = tc.addStructure(name=\"dummy\", a=3, b=2, c=1.5, alpha=50, beta=70, gamma=60)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vec = tc.miller(0, 0, 1)\n",
    "xaxis = np.array([1, 0, 0])\n",
    "listvec2 = np.array([[1, 2, 3],[4, 5, 6]]).T\n",
    "listvec3 = np.array([[1, 2, 3],[4, 5, 6], [7, 8, 9]]).T\n",
    "steel.realToRecyp(listvec2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "jeol3000 = tc.addMicroscope(kv=300)\n",
    "doubletilt = jeol3000.addStage(name=\"doubletilt\")\n",
    "ccd = jeol3000.addDetector(name=\"MSC\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ccd.setCalibration(\"rotationCalibrationMSC-IMG.txt\", mode = \"imaging\", type = \"r\")\n",
    "ccd.setCalibration(\"rotationCalibrationMSC-DIFF.txt\", mode = \"diffraction\", type = \"r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "doubletilt.setAlpha(0)\n",
    "doubletilt.setBeta(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#ccd.plotRealAxes(\"i\", 150000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "crl = doubletilt.addCrystal(\"steel\", name = \"faggot\")\n",
    "crl2 = doubletilt.addCrystal(\"dummy\", name = \"dickins\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "crl.calcOrient(tc.miller(1, 1, 0), tc.miller(1, -1, 1), -30, \"MSC\", \"d\", 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "crl.plotCrystalonDetector(\"MSC\", \"d\", 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "testcase = np.array([[0,2,16],[1,0,8],[-1, -1, 8]]).T\n",
    "#test if calcAlphaBeta works\n",
    "doubletilt.calcAlphaBeta(testcase)\n",
    "#it works. Note that the the vector needs to be rotated to the vector INTO the plane!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#now lets change the tilt of the stage\n",
    "doubletilt.setAlpha(10)\n",
    "doubletilt.setBeta(20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#note that the points plotted are the ones pointing DOWN. The zone is \n",
    "#whatever is pointing into the page\n",
    "crl.plotCrystalonDetector(\"MSC\", \"d\", 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#check if the conversions are working correctly\n",
    "\n",
    "chk = tc.miller(1, 1, 0)\n",
    "print(crl.millerToDetector(chk, \"MSC\", \"d\", 30))\n",
    "print(crl.detectorToMiller(xaxis, \"MSC\", \"d\", 30))\n",
    "print(crl.detectorToMiller(tc.sizeAngle(1, 30), \"MSC\", \"d\", 30, typ=\"recyp\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#crl.millerToDetector(listvec2, \"MSC\", \"i\", 60000)\n",
    "#zax = crl.detectorToMiller(tc.miller(0,0,3.66), \"MSC\", \"i\", 60000, verbose = True)\n",
    "#xax = crl.detectorToMiller(tc.miller(np.cos(30/360*2*np.pi),np.sin(30/360*2*np.pi),0), \"MSC\", \"i\", 60000, verbose = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#see how the tilting is working\n",
    "crl.getAlphaBetaMiller(tc.miller(1,1,0), verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#This seems to work\n",
    "#what if we first have some angle, then set the orientation, can we still calculate it?\n",
    "#now lets change the tilt of the stage\n",
    "doubletilt.setAlpha(10)\n",
    "doubletilt.setBeta(20)\n",
    "\n",
    "crl.calcOrient(tc.miller(1, 1, 0), tc.miller(1, -1, 1), -30, \"MSC\", \"d\", 30)\n",
    "\n",
    "crl.plotCrystalonDetector(\"MSC\", \"d\", 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "crl.getAlphaBetaMiller(tc.miller(1,1,0), verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#it works! What about all combinations?\n",
    "all110 = crl.getSym(tc.miller(1, 1, 0), checksame=True, typ = \"real\")\n",
    "print(all110)\n",
    "crl.getAlphaBetaMiller(all110)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#let's change the range so we see more posibilities\n",
    "doubletilt.setAlphaRange(-90, 90)\n",
    "doubletilt.setBetaRange(-90, 90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "crl.getAlphaBetaMiller(all110).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#test it out\n",
    "doubletilt.tiltTo((2.01, -39.87))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "crl.plotCrystalonDetector(\"MSC\", \"i\", 4e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "crl.getZone(verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "steel.getVectorLength(tc.miller(2, 0, 0), typ=\"recyp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "##Check if the two-beam tilting works\n",
    "crl.calcSampleTilt(tc.miller(2, 0, 0), n=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "jeol3000.setKv(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "##check if the two-beam tilting works if we change the accelerating voltage\n",
    "crl.calcSampleTilt(tc.miller(2, 0, 0), n=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "def integerRep(vec, decimals = 9):\n",
    "    \"\"\"Make a vector or array of column vectors more presentable by giving their closest integer representation.\"\"\"\n",
    "    vec = np.array(vec)\n",
    "    vec = np.round(vec, decimals = decimals)\n",
    "    \n",
    "    def lcm(a, b):\n",
    "        a = int(a)\n",
    "        b = int(b)\n",
    "        return a*b/math.gcd(a,b)\n",
    "\n",
    "    def lcm3(lst):\n",
    "        return lcm(lst[0], lcm(lst[1], lst[2]))\n",
    "\n",
    "    def getDenominator(x):\n",
    "        return Fraction(str(x)).limit_denominator().denominator\n",
    "    \n",
    "    def getIntRep(vect):\n",
    "        print(vect)\n",
    "        vect = vect/np.min(np.abs(vect[np.nonzero(vect)]))\n",
    "        g=[]\n",
    "        for i in vect.tolist():\n",
    "            i = abs(i)\n",
    "            if i!=0.0:\n",
    "                g.append(getDenominator(i))\n",
    "            else:\n",
    "                g.append(1)\n",
    "        return vect*lcm3(g)\n",
    "    \n",
    "    if vec.ndim ==1:\n",
    "        return getIntRep(vec)\n",
    "        \n",
    "    elif vec.ndim ==2:\n",
    "        toret = np.zeros(vec.shape)\n",
    "        for i in range(vec.shape[1]):\n",
    "            vecy = vec[:,i]\n",
    "            toret[:,i] = getIntRep(vecy)\n",
    "        return toret   \n",
    "    \n",
    "    else:\n",
    "        print(\"The data is not supplied in the right format\")\n",
    "        return vec\n",
    "\n",
    "#Non working method to calculate the intersection with the reflection axis with the ewald sphere\n",
    "def calcNfromAngles(g, sa, ba, K0 = 507.93):\n",
    "    \"\"\"This method calculates the active 'n' after a given sample angle tilt and beam angle tilt\"\"\"\n",
    "    #get the length of vector g in nm-1\n",
    "    length = g\n",
    "    #get the length of the radius of the ewald sphere in nm-1\n",
    "    #K0 = 507.93\n",
    "    #takes sa and ba as radian quantitites!\n",
    "    #this calculation was done on paper, sa = sample angle = phi, ba = beam angle = theta\n",
    "    c = np.cos(ba) - np.sin(ba)*np.tan(sa)\n",
    "    a = np.tan(sa)\n",
    "    b = -1\n",
    "    #solve the equation a*cos(v)+b*sin(v) = c\n",
    "    v = np.arccos(c/K0) + np.arccos(a/K0)\n",
    "    #now get x and y\n",
    "    x = K0*np.cos(v) + K0*np.cos(ba)\n",
    "    y = x*np.tan(sa)\n",
    "    #now get the length of this vector\n",
    "    l = np.sqrt(x**2+y**2)\n",
    "    #now divide by the length of g to get n\n",
    "    return l/length\n",
    "    \n",
    "#implemented in the file\n",
    "def sampleTilt(g, n, units = \"degrees\", K0 = 507.93):\n",
    "    \"\"\"Sample tilt necessary to put reflection ng into two-beam condition\"\"\"\n",
    "    #get the length of vector g in nm-1\n",
    "    length = g\n",
    "    #get the length of the radius of the ewald sphere\n",
    "    #convert to nm-1\n",
    "    #K0 = 507.93\n",
    "    ans = np.pi/2 - np.arccos(n*length/2/K0)\n",
    "    if units == \"degrees\":\n",
    "        return ans/(2*np.pi)*360\n",
    "    elif units == \"mrad\":\n",
    "        return ans*1000\n",
    "    elif units == \"radians\":\n",
    "        return ans\n",
    "    else:\n",
    "        print(\"Units not recognized\")\n",
    "        return 0\n",
    "        \n",
    "\n",
    "    \n",
    "def beamTilt(g, n, k, units = \"degrees\", K0 = 507.93):\n",
    "    \"\"\"This calculates how much the beam needs to be tilted when beam ng is active in order for beam kg to get on the optical axis\"\"\"\n",
    "    #get the length of vector g in nm-1\n",
    "    length = g\n",
    "    #get the length of the radius of the ewald sphere\n",
    "    #convert to nm-1\n",
    "    #K0 = 507.93\n",
    "    yot = sampleTilt(g, n, units = \"radians\", K0=K0)\n",
    "    thet = np.pi/2 - np.arccos(k*length*np.cos(yot)/K0)\n",
    "    if units == \"degrees\":\n",
    "        return thet/(2*np.pi)*360\n",
    "    elif units == \"mrad\":\n",
    "        return thet*1000\n",
    "    elif units == \"radians\":\n",
    "        return thet\n",
    "    else:\n",
    "        print(\"Units not recognized\")\n",
    "        return 0\n",
    "    \n",
    "def calcNfromAngles(g, sa, ba, K0 = 507.93):\n",
    "    \"\"\"This method calculates the active 'n' after a given sample angle tilt and beam angle tilt\"\"\"\n",
    "    #get the length of vector g in nm-1\n",
    "    length = g\n",
    "    #get the length of the radius of the ewald sphere in nm-1\n",
    "    #K0 = 507.93\n",
    "    #takes sa and ba as radian quantitites!\n",
    "    #this calculation was done on paper, sa = sample angle = phi, ba = beam angle = theta. This was purely geometric.\n",
    "    ang = np.pi/2-sa-ba\n",
    "    dst = 2*K0*np.cos(ang)\n",
    "    return dst/length\n",
    "\n",
    "def calcNfromTilts(g, n, k, K0 = 507.93):\n",
    "    \"\"\"This method calculates what the active 'n' value is after a sample tilt and beamtilt\"\"\"\n",
    "    yot = sampleTilt(g, n, units = \"radians\", K0=K0)\n",
    "    thet = beamTilt(g, n, k, units = \"radians\", K0 = K0)\n",
    "    ans = calcNfromAngles(g, yot, thet, K0 = K0)\n",
    "    return ans\n",
    "\n",
    "\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "def calcWB2beam(g, k, i, K0 = 507.93):\n",
    "    \"\"\"This method calculates what the correct 2-beam tilt condition should be (n) if we want the ewald sphere to cut through ig after the beam tilt so that reflection k is centered\"\"\"\n",
    "    #get the length of vector g\n",
    "    length = g\n",
    "    #the angle between the beam and the reflections can be easily calculated \n",
    "    ang = np.arccos(i*length/2/K0)\n",
    "    #need to solve the following equation numericaly:\n",
    "    def eq(phi):\n",
    "        return np.arccos(k*length*np.cos(phi)/K0)-ang-phi\n",
    "    \n",
    "    phi = fsolve(eq, 0)\n",
    "    #from the known phi, the initial 2-beam can be determined\n",
    "    lent = 2*K0*np.sin(phi)\n",
    "    \n",
    "    return (lent/length)[0]\n",
    "\n",
    "def calcSg(g, n, k, n2, K0 = 507.93, verbose=False):\n",
    "    #n and k determine the setup conditions: sample tilt and beam tilt\n",
    "    #get the length of vector g\n",
    "    length = g\n",
    "    #n2 is used for calculating sg, usually this will be k, the reflection used for imaging \n",
    "    yot = sampleTilt(g, n, units = \"radians\", K0=K0)\n",
    "    thet = beamTilt(g, n, k, units = \"radians\", K0 = K0)\n",
    "    #slope and intercept of the line\n",
    "    slop = np.tan((np.pi-thet)/2)\n",
    "    intc = n2*length*(np.sin(yot)-slop*np.cos(yot))\n",
    "    #solve the intersection with the circle\n",
    "    c = intc/K0 + slop*np.sin(thet)-np.cos(thet)\n",
    "    a = -slop\n",
    "    b = 1\n",
    "    #v = np.arcsin(c/np.sqrt(a**2+b**2))-np.arcsin(a/np.sqrt(a**2+b**2))\n",
    "    fp = np.arcsin(-c/np.sqrt(a**2+b**2))\n",
    "    sp = np.arcsin(a/np.sqrt(a**2+b**2))\n",
    "    v1 = fp - sp - np.pi\n",
    "    v2 = -fp - sp\n",
    "    \n",
    "    if np.sin(v1)<0:\n",
    "        v = v1\n",
    "    elif np.sin(v2)<0:\n",
    "        v = v2\n",
    "    else:\n",
    "        print(\"Something went terribly wrong.\")\n",
    "        v=0\n",
    "    \n",
    "    #print(fp/(2*np.pi)*360)\n",
    "    #print(sp/(2*np.pi)*360)\n",
    "    #print(ad/(2*np.pi)*360)\n",
    "    \n",
    "    #what are the points on the ewald sphere\n",
    "    x = K0*(np.cos(v)+np.sin(thet))\n",
    "    y = K0*(np.sin(v)+np.cos(thet))\n",
    "    #choose the point where y is lower than the center of the ewald sphere\n",
    "    #location of x and y of n2g\n",
    "    x0 = n2*length*np.cos(yot)\n",
    "    y0 = n2*length*np.sin(yot)\n",
    "    \n",
    "    if verbose:\n",
    "        print(\"v: %s\" %(v/np.pi/2*360))\n",
    "        print(\"x_gn2: %s\" %x0)\n",
    "        print(\"y_gn2: %s\" %y0)\n",
    "        print(\"slope: %s\" %slop)\n",
    "        print(\"intercept: %s\" %intc)\n",
    "        print(\"x_int: %s\" %x)\n",
    "        print(\"y_int: %s\" %y)\n",
    "    \n",
    "    #the distance:\n",
    "    if y0<y: #the ewald sphere is above the n2g point - the point is outside - negative sg\n",
    "        return -np.sqrt((x0-x)**2+(y0-y)**2)\n",
    "    else:\n",
    "        return np.sqrt((x0-x)**2+(y0-y)**2)\n",
    "\n",
    "\n",
    "####Non working way to solve for the tilts of the holder\n",
    "\n",
    "from scipy.optimize import fsolve\n",
    "import math\n",
    "\n",
    "toconv = np.array([0, 1, 0])\n",
    "\n",
    "def equations(p):\n",
    "    alpha, beta, dummy = p\n",
    "    c = toconv[0]\n",
    "    d = toconv[1]\n",
    "    f = toconv[2]\n",
    "    L = np.sqrt(c**2+d**2+f**2)\n",
    "    \n",
    "    eq1 = c*np.cos(beta) + f*np.sin(beta)\n",
    "    eq2 = c*np.sin(alpha)*np.sin(beta) - f*np.sin(alpha)*np.cos(beta) + d*np.cos(alpha)\n",
    "    eq3 = -c*np.sin(beta)*np.cos(alpha) + f*np.cos(alpha)*np.cos(beta) + d*np.sin(alpha)-L\n",
    "    return eq1, eq2, eq3\n",
    "\n",
    "\n",
    "\n",
    "alp, bet, dummy =  fsolve(equations, (1, 1, 1))\n",
    "\n",
    "\"\"\"\n",
    "if abs(abs(alp)-np.pi)<1e-10:\n",
    "    alp = 0\n",
    "if abs(abs(bet)-np.pi)<1e-10:\n",
    "    bet = 0\n",
    "\"\"\"\n",
    "    \n",
    "print((alp/(2*np.pi)*360, bet/(2*np.pi)*360, dummy))\n",
    "print((round(alp/np.pi/2*360,5)%(360), round(bet/np.pi/2*360,5)%(360), dummy))\n",
    "print(equations((alp%(2*np.pi), bet%(2*np.pi), dummy)))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "tester = np.array([[-1, -1, 8],[2, 3, 4],[4, 2, 6],[1, 3, 5]]).T\n",
    "\n",
    "def sak(vec):\n",
    "    return np.array([vec[0] + vec[2], vec[0] + vec[1]]) \n",
    "\n",
    "def sak2(vec):\n",
    "    return vec[0] + vec[2]\n",
    "\n",
    "def vecto(func, lstvc, **kwargs):\n",
    "    cols = np.size(lstvc, 1)\n",
    "    ans = np.array([])\n",
    "    for i in range(cols):\n",
    "        act = lstvc[:,i]\n",
    "        #add the first entry\n",
    "        if not ans.size:\n",
    "            ans = np.array(func(act, **kwargs))\n",
    "        else:\n",
    "            ans = np.vstack([ans, np.array(func(act, **kwargs))])\n",
    "    #if the answer is only one column, return a one dimensional array\n",
    "    if np.size(ans, 1)==1:\n",
    "        return ans[:,0]\n",
    "    else: #otherwise return a matrix\n",
    "        return ans.T"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "from fractions import Fraction\n",
    "import math\n",
    "\n",
    "tester = tc.miller(123, 43, 11)\n",
    "#testn = np.array(tc.normMatrix(tester))\n",
    "testn = testn/np.min(testn)\n",
    "#testn = testn[:, 0]\n",
    "print(testn)\n",
    "\n",
    "def lcm(a, b):\n",
    "    a = int(a)\n",
    "    b = int(b)\n",
    "    return a*b/math.gcd(a,b)\n",
    "\n",
    "def lcm3(lst):\n",
    "    return lcm(lst[0], lcm(lst[1], lst[2]))\n",
    "\n",
    "def getDenominator(x):\n",
    "    return Fraction(str(x)).limit_denominator().denominator\n",
    "\n",
    "g=[]\n",
    "for i in testn.tolist():\n",
    "    g.append(getDenominator(i))\n",
    "    \n",
    "testn*lcm3(g)\n",
    "\n",
    "def inrange(a,b):\n",
    "    a=a/(2*np.pi)*360\n",
    "    b=b/(2*np.pi)*360\n",
    "    atop = 30\n",
    "    alow = -30\n",
    "    btop = 25\n",
    "    blow = -25\n",
    "    if a<=atop and a>=alow and b<=btop and b>=blow:\n",
    "        return True\n",
    "    else:\n",
    "        return False"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "d001 = crl.getSym(tc.miller(1, 0, 0), typ=\"real\")\n",
    "toplotd001 = crl.millerToDetector(d001, \"MSC\", \"i\", 100000, typ=\"real\", verbose=True)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "tc.normMatrix(toplotd001)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "test = tc.miller(1.5, 2.3, 3.2)\n",
    "test = crl.getSym(test)\n",
    "\n",
    "testlab = tc.miller(1, 3, 6)\n",
    "testlab = crl.getSym(testlab)\n",
    "\n",
    "fig, ax = plt.subplots(1)\n",
    "tc.stereographicCanvas(ax)\n",
    "tc.plotStereographic(test, testlab, ax)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "def continuedFraction(lst):\n",
    "    a=1/lst[-1]\n",
    "    inx = 1\n",
    "    while(inx<len(lst)-1):\n",
    "        print(lst[-1-inx])\n",
    "        a=1/(lst[-1-inx]+a)\n",
    "        inx+=1\n",
    "    a = a+lst[0]\n",
    "    return a\n",
    "\n",
    "continuedFraction(np.array([1, 1, 1, 3]))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "def calcForOne(toz, rnd = 10, verbose=False):\n",
    "    c = toz[0]#+1e-15\n",
    "    d = toz[1]#+1e-15\n",
    "    f = toz[2]#+1e-15\n",
    "\n",
    "    length = np.sqrt(c**2+d**2+f**2)\n",
    "\n",
    "    #there is a 180 degree uncertainty on the answer. If the angle is negative (fourth quadrant) it could be meant for the second quadrant.\n",
    "    beta = np.arctan(-c/f)\n",
    "    beta2 = 0\n",
    "    if beta<0:\n",
    "        beta2 = np.pi + beta\n",
    "    else:\n",
    "        beta2 = -np.pi + beta\n",
    "    \n",
    "    #there is an exception: for all vectors in the y-z plane (c=0), beta can be any value. So it is best to set it to 0\n",
    "    if round(c,6)==0:\n",
    "        beta=0\n",
    "        beta2=np.pi\n",
    "    \n",
    "    #this results in two possibilities for the solution of the second equation\n",
    "    denom1 = f*np.cos(beta) - c*np.sin(beta)\n",
    "    denom2 = f*np.cos(beta2) - c*np.sin(beta2)\n",
    "\n",
    "    alpha1a = np.arctan(d/denom1)\n",
    "    alpha1b = 0\n",
    "    if alpha1a<0:\n",
    "        alpha1b = np.pi + alpha1a\n",
    "    else:\n",
    "        alpha1b = -np.pi + alpha1a\n",
    "\n",
    "    alpha2a = np.arctan(d/denom2)\n",
    "    alpha2b = 0\n",
    "    if alpha2a<0:\n",
    "        alpha2b = np.pi + alpha2a\n",
    "    else:\n",
    "        alpha2b = -np.pi + alpha2a\n",
    "\n",
    "    #the third equation serves as a test - this needs to be 0 for a correct solution\n",
    "    def test(a, b):\n",
    "        return (np.round(-c*np.sin(b)*np.cos(a) + f*np.cos(a)*np.cos(b) + d*np.sin(a) - length, 5) == 0 and inrange(a,b))\n",
    "\n",
    "    possible = []\n",
    "    #test all the possibilities\n",
    "    if test(alpha1a, beta):\n",
    "        if verbose:\n",
    "            print(\"Option 1 is possible\")\n",
    "        possible.append([np.round(alpha1a/(2*np.pi)*360, rnd), np.round(beta/(2*np.pi)*360, rnd)])\n",
    "    if test(alpha1b, beta):\n",
    "        if verbose:\n",
    "            print(\"Option 2 is possible\")\n",
    "        possible.append([round(alpha1b/(2*np.pi)*360, rnd), round(beta/(2*np.pi)*360, rnd)])\n",
    "    if test(alpha2a, beta2):\n",
    "        if verbose:\n",
    "            print(\"Option 3 is possible\")\n",
    "        possible.append([round(alpha2a/(2*np.pi)*360, rnd), round(beta2/(2*np.pi)*360, rnd)])\n",
    "    if test(alpha2b, beta2):\n",
    "        if verbose:\n",
    "            print(\"Option 4 is possible\")\n",
    "        possible.append([round(alpha2b/(2*np.pi)*360, rnd), round(beta2/(2*np.pi)*360, rnd)])\n",
    "\n",
    "    return possible\n",
    "\n",
    "tester = np.array([-1, -1, 8])\n",
    "calcforone(tester, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}