#!/usr/bin/env python
# -*- coding: utf-8 -*-

from numpy import *
#from matplotlib.pyplot import *
from numpy.linalg import norm
from scipy.special import ellipk, ellipe
#import time


#  functions for the calculation of the magnetic field around the current loop, by function Bloop_integ 
#is values calculated by integration, in the function Bloop_analytic is used analytical formula
#Autor :Tomas Odstrcil


# obecný pruodový profil může být nahrazen skupinou bodých proudových profilů  v těžištích příslušných mezikruží nesoucích proud 

global mu_0,a

a = 0.085		#[m]
R_0 = 0.4		#[m]
mu_0 = 4*pi*1e-7
p = 2		#píkovací faktor


    
def CurrentProfile(rho,zeta,Rc,zc,I, Shaf_shift, centered = 'False'):		
#parabolický profil, rho,zeta jsou souřadnice, Rc,zeta - střed plazmatu (ale ne proudového profilu!),
#Shaf_shift veličina udávají shafranův posuv  většinou je 0-0.2
    j= zeros(shape(rho)+shape(Rc)) #předpokládám že jeden z nich je číslo
    DRc = Rc-R_0
    Drho = rho-R_0
    ao = a-hypot(DRc,zc)
    n = size(rho,0)

    if ao/a <=5./n:  #plasma escaped out of the chamber
	DRc/= hypot(DRc,zc)*(1-5./n)*a
	zc/= hypot(DRc,zc)*(1-5./n)*a
	ao = 5*a/n

    if centered:
	r2 = (zeta)**2+(Drho)**2   #centered profile
    else:
	r2 = (zeta-zc)**2+(Drho-DRc)**2   #shifted profile
    plazma = (r2 < ao**2)

    j[plazma] = (p+1)*I/(pi*ao**2)*(1-r2[plazma]/ao**2)**p
    j[plazma]*= (1+(j[plazma])/amax(j[plazma])*Shaf_shift*sin(pi*(rho[plazma]-Rc)/ao))

    return j,ao
    
    
    
def Bloop_integ(R,z,rho,zeta):	# toroidální vodič - vypočtenu přímou integrací
    #R,Z poloha detektoru, rho, zeta poloha vlákna

    m = 1000
    dphi = pi/m  #protože jsem integroval jen přes půlkruh
    phi_i = (arange(0,m))*dphi
    B = zeros(3)
    R_Q = zeros(3)
    R_Q[2] = zeta
    dR_Q = zeros(3)
    dR = zeros(3)
    R_D = zeros(3)
    R_D[0] = R
    R_D[2] = z
    S = sin(phi_i)
    C = cos(phi_i)
    N = arange(m)
    Bvec= empty((m,3))
    for i,c,s in zip(N,C,S):

	
	R_Q[0] = rho*c
	R_Q[1] = rho*s
	dR_Q[0] = -s
	dR_Q[1] = c
	dR = R_D-R_Q
	dB = cross(dR_Q, dR)/norm(dR)**3
	Bvec[i,:] = dB
	if i == 0 or i == m:   #simpson integration formula
	    B+= dB
	elif i%2 == 0:
	    B+= 2*dB
	else:
	    B+= 4*dB
    
     
    B*=mu_0/(4*pi)*rho*(dphi/3)
    B*= 2  #protože jsme integroval jen přes půlkruh
    B[1] = 0#protože jsme integroval jen přes půlkruh
    Br = B[0]
    Bz = B[2]
    return Bz, Br      
    

def Bloop_analytic(R,z,rho,zeta,out=None):	
    #R,Z poloha detektoru, rho, zeta poloha vlákna

    #http://www.netdenizen.com/emagnettest/offaxis/?offaxisloop

    r = abs(R)
    x = z-zeta

    
    B0 = mu_0/(2.*pi*rho)
    alpha = divide(r,rho)
    gamma = divide(x,r,r)
    beta_sq = square(divide(x,rho,x),x)
    Q = square(alpha+1.)+beta_sq
    sQ = sqrt(Q)
    k_sq = divide(alpha,Q)*4.
    
    

    E_val =  ellipe(k_sq)
    K_val =  ellipk(k_sq)
    
    
    if out == None:
	B = empty(shape(rho)+shape(R)+ (2,)) # předpokládám že jedno z nich  je pole a druhé konstanta
    else:
	B = out

    Bz = (E_val*(1.-square(alpha)-beta_sq)/(Q-alpha*4.)+K_val)*B0/sQ# B_z
    Br = B0*gamma/sQ*(E_val*(Q-alpha*2.)/(Q-alpha*4.)-K_val)

    B[...,1] = Bz
    B[...,0] = Br

    return B