gravcaustic.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 16 19:50:54 2021

caustic.py

@author: nolte

D. D. Nolte, Optical Interferometry for Biology and Medicine (Springer,2011)
"""

import numpy as np
from matplotlib import pyplot as plt
from numpy import random as rnd
from scipy import signal as signal

plt.close('all')

N = 256

def grav(sy,sx,w):
    
    x = np.arange(-sx/2,sy/2,1)
    y = np.arange(-sy/2,sy/2,1)
    y = y[..., None]
    
    ex = np.ones(shape=(sy,1))
    x2 = np.kron(ex,x**2);
    
    ey = np.ones(shape=(1,sx));
    y2 = np.kron(y**2,ey);

    rad2 = (x2+y2);

    A = np.positive(1.0/np.sqrt((1.0 - w/np.sqrt(rad2)))**2 + 1e-10) ;

    return A

def gauss2(sy,sx,wy,wx):
    
    x = np.arange(-sx/2,sy/2,1)
    y = np.arange(-sy/2,sy/2,1)
    y = y[..., None]
    
    ex = np.ones(shape=(sy,1))
    x2 = np.kron(ex,x**2/(2*wx**2));
    
    ey = np.ones(shape=(1,sx));
    y2 = np.kron(y**2/(2*wy**2),ey);

    rad2 = (x2+y2);

    A = np.exp(-rad2);

    return A

def cauchy(sy,sx,wy,wx):
     
    x = np.arange(-sx/2,sy/2,1)
    y = np.arange(-sy/2,sy/2,1)
    y = y[..., None]
    
    ex = np.ones(shape=(sy,1))
    x2 = np.kron(ex,x**2/(2*wx**2))
    
    ey = np.ones(shape=(1,sx))
    y2 = np.kron(y**2/(2*wy**2),ey)

    rad2 = (x2+y2)

    A = 1/(1+rad2)

    return A

def speckle2(sy,sx,wy,wx):

    Btemp = 2*np.pi*rnd.rand(sy,sx);
    B = np.exp(complex(0,1)*Btemp);

    C = gauss2(sy,sx,wy,wx);

    Atemp = signal.convolve2d(B,C,'same');

    Intens = np.mean(np.mean(np.abs(Atemp)**2));

    D = np.real(Atemp/np.sqrt(Intens));

    Dphs = np.arctan2(np.imag(D),np.real(D));

    return D, Dphs


#Sp = grav(N,N,N/8)

Sp = cauchy(N,N,N/8,N/8)

plt.figure(2)
plt.matshow(Sp,2,cmap=plt.cm.get_cmap('seismic'))  # hsv, seismic, bwr
plt.show()

fx, fy = np.gradient(Sp);

fxx,fxy = np.gradient(fx);
fyx,fyy = np.gradient(fy);

J = fxx*fyy - fxy*fyx;

D = np.abs(1/J)

plt.figure(3)
plt.matshow(D,3,cmap=plt.cm.get_cmap('gray'))  # hsv, seismic, bwr
plt.clim(0,0.5e7)
plt.show()

eps = 1e-7
cnt = 0
E = np.zeros(shape=(N,N))
for yloop in range(0,N-1):
    for xloop in range(0,N-1):
        
        d = 5*N/2
        
        indx = int(N/2 + (d*(fx[yloop,xloop])+(xloop-N/2)/2))
        indy = int(N/2 + (d*(fy[yloop,xloop])+(yloop-N/2)/2))
        
        if ((indx > 0) and (indx < N)) and ((indy > 0) and (indy < N)):
            E[indy,indx] = E[indy,indx] + 1

plt.figure(4)
plt.imshow(E,interpolation='bicubic',cmap=plt.cm.get_cmap('gray'))
plt.clim(0,20)
plt.xlim(N/4, 3*N/4)
plt.ylim(N/4,3*N/4)