week3.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "week3.ipynb",
      "provenance": [],
      "authorship_tag": "ABX9TyMmxT2fNpy/PJ9THulhCySx",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/xspronken/bigdataW1/blob/master/week3.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oji3R7Zokxoy",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import math as m\n",
        "\n",
        "\n",
        "''' For conversion from natural units'''\n",
        "G = 6.67*10**(-11)\n",
        "c = 3*10**8\n",
        "\n",
        "''' defining Power Spectral Density S(f) '''\n",
        "fs = 20 #lower limit\n",
        "f0 = 215 #upper limit\n",
        "\n",
        "def S(x):\n",
        "  s0 = 10**(-49)\n",
        "  return s0*((4.49*x)**(-56)+0.16*x**(-4.52)+0.52+0.32*x**2)\n",
        "\n",
        "'''variables'''\n",
        "SNR = 10\n",
        "def A2(SNR, dx):\n",
        "  accumulator = 0\n",
        "  for i in range (int(fs/dx), int(f0/dx)):\n",
        "   accumulator += (i*dx)**(-7/6) / (S(i*dx))\n",
        "  return SNR/(4*accumulator)\n",
        "A = m.sqrt(A2(SNR, 0.001))\n",
        "\n",
        "m1 = 10**30/((c**2)/G)\n",
        "m2 = 10**30/((c**2)/G)\n",
        "f = 150\n",
        "\n",
        "M = (m1+m2)\n",
        "n = m1*m2/(M**2)\n",
        "t_c = 100\n",
        "phi_c = 100\n",
        "\n",
        "params = [t_c, phi_c, M, n]\n",
        "\n",
        "alpha = [1, 0, (20/9*(743/336 + n*11/4)), -16*m.pi, 10*(3058673/1016064 + n*5429/1008 + (n**2) * 617/144)]\n",
        "PN = 1\n",
        "\n",
        "\n",
        "''' defining the Fourrier transform h(f)'''\n",
        "def h(f, params, x):\n",
        "  #defining the waveform's variables\n",
        "  t_c, phi_c, M , n = params\n",
        "  v = (m.pi*M*f)**(1/3)\n",
        "  \n",
        "  #sum of all coefficients (alpha) in the fourrier phase\n",
        "  x *= 2\n",
        "  sum_alpha = 0\n",
        "  for i in range(x):\n",
        "    sum_alpha += (alpha[i] * (v**i))\n",
        "  \n",
        "  #defining the waveform psi(f)\n",
        "  psi =  2 * m.pi * f * t_c - phi_c - m.pi/4 + 3/(128*n*v**5) * sum_alpha\n",
        "  \n",
        "  return  A * f**(-7/6) * m.e**(1j * psi)\n",
        "\n",
        "\n",
        "''' Calculating the partial derivative of h(f) '''\n",
        "def h_partial(f, params, x, dx):   \n",
        "  dparams = params.copy()\n",
        "  dparams[x] += dx\n",
        "  dy = h(f, dparams, PN)-h(f, params, PN)\n",
        "  return dy/dx\n",
        "\n",
        "''' Generating the Fischer information matrix '''\n",
        "def generate_matrix(params, dx):\n",
        "  \n",
        "# create an n*n matrix where n is the number of observables\n",
        "  n = len(params)\n",
        "  matrix = np.empty([n,n]) \n",
        "  \n",
        "# Start filling the matrix\n",
        "  for i in range (n):\n",
        "    for j in range (n):\n",
        "      \n",
        "      #start integration\n",
        "      accumulator = 0\n",
        "      for k in range(int(fs/dx), int(f0/dx)):\n",
        "        part_i = h_partial(k*dx, params,i, 0.001)\n",
        "        part_j = np.conj(h_partial(k*dx, params,j, 0.001))\n",
        "        accumulator += np.real(part_i * part_j / S(k*dx))\n",
        "      \n",
        "      #once integration is done add it to the matrix\n",
        "      matrix[i][j] = accumulator*4\n",
        "  \n",
        "  return matrix\n",
        "\n",
        "\n",
        "''' Extract the diagonal elements from the inverted Fischer matrix'''\n",
        "\n",
        "def extract_from_matrix(matrix):\n",
        "  #invert matrix\n",
        "  y = np.linalg.inv(matrix) \n",
        "  z = []\n",
        "  \n",
        "  #fetch the square root of diagonal elements\n",
        "  for i in range (len(y)):\n",
        "    z.append(m.sqrt(y[i][i]))\n",
        "  \n",
        "  return z\n",
        "\n",
        "test = extract_from_matrix(generate_matrix(params, 0.001))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gZKgeXY_Y0UV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "test"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}