# -*- coding: utf-8 -*- import numpy as np import cmath as cm from math import * from decimal import * ############################################################################################# ############################################################################################# ########## just call the potential with ##################################################### ##### def negf(xx): # inverse of the fitness aka negf is the potential -- to minimize ####### ###### you can have upto 5 variables in total, but easily extendable ######################## ############################################################################################# ############################################################################################# ###### The following are to make it easy to copy the function from Mathematica using CForm ## ######## It's a pain that it gives lines that end in / <-- you have to remove the new line ## def Power(c,s): # function to enable copying from mathematica using CForm[func] return np.power(c,s) def Sqrt(a): # function to enable copying from mathematica using CForm[func] return np.power(a,0.5) def Cos(c): # function to enable copying from mathematica using CForm[func] return np.cos(c) def Sin(c): # function to enable copying from mathematica using CForm[func] return np.sin(c) ############## RESPOSITORY OF FUNCTIONS ########## ############## IN THE FORM ######################## # names_myname = [r'$\tau_s$', r'$\log_{10}(\tau_b)$', r'$\log_{10}(vol)$', r'$\rho_1$', r'$\rho_2$'] # def negf_myname(xx): ################################################### ############## names_wobbly_gaussian = [r'$x$',r'$y$', r'dummy', r'dummy', r'dummy'] def negf_wobbly_guassian(xx): # inverse of the fitness aka the potential -- to minimize x0 = np.array(xx[0]) x1 = np.array(xx[1]) # V = -(36*np.sin(2*x1)*np.cos(2*x0) + 12*(x0 + x1) - Power(x1,2) - Power(x0,2) ) V = -(np.exp(-0.1*(x0**2+x1**2)) * np.sin(2*x1)*np.cos(2*x0) ) return V ############## names_racetrack = [r'$\tau_1$',r'$\tau_2$', r'$a_1\rho_1/\pi $', r'$a_2\rho_2/\pi $', r'dummy'] def negf_racetrack(params,xx): # inverse of the fitness aka the potential -- to minimize (here treating as a 2d system with set angles global content2 E=2.7182818284590452354 Pi=3.1415926535897932385 t1 = np.array(xx[0]) t2 = np.array(xx[1]) r1 = 0 r2 = 0 Kcs=params[0] k111=params[1] chiiP11169=params[2] W0=params[3] A1=params[4] a1=params[5] A2 =params[6] a2=params[7] gs=params[8] xi = -1.2*chiiP11169/(16.0 * pi**3) xihat = xi/Power(gs,3.0/2.0) Vc = ( (2*Power(E,Kcs)*gs*((3*Power(W0,2)*xihat*(9*Power(t1,6) - 18*Power(t1,5)*t2 + 15*Power(t1,4)*Power(t2,2) - 6*Power(t1,3)*(Power(t2,3) - 7*xihat) + 14*t1*Power(t2,2)*xihat + 4*Power(xihat,2) + Power(t1,2)*(Power(t2,4) - 42*t2*xihat)))/ (2.*(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) - 2*xihat)* Power(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat,2)) + Power(A2,2)*Power(E,a2*t1*(3*t1 - 2*t2))* ((Power(a2,2)*Power(t1,2)*Power(3*t1 - 2*t2,2)* (6*Power(t1,3) - 6*Power(t1,2)*t2 + 2*t1*Power(t2,2) - xihat))/ (6*Power(t1,3) - 6*Power(t1,2)*t2 + 2*t1*Power(t2,2) - 4*xihat) - 2*Power(a2,2)*t1*(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat) + (3*xihat*(9*Power(t1,6) - 18*Power(t1,5)*t2 + 15*Power(t1,4)*Power(t2,2) - 6*Power(t1,3)*(Power(t2,3) - 7*xihat) + 14*t1*Power(t2,2)*xihat + 4*Power(xihat,2) + Power(t1,2)*(Power(t2,4) - 42*t2*xihat)))/ (2.*(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) - 2*xihat)* Power(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat,2)) - (a2*t1*(3*t1 - 2*t2)*(18*Power(t1,6) - 36*Power(t1,5)*t2 + 30*Power(t1,4)*Power(t2,2) + t1*Power(t2,2)*xihat + 8*Power(xihat,2) + 3*Power(t1,3)*(-4*Power(t2,3) + xihat) + Power(t1,2)*(2*Power(t2,4) - 3*t2*xihat)))/ ((3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) - 2*xihat)* (3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat))) + (Power(A1,2)*((Power(a1,2)*Power(-3*t1 + t2,4)* (6*Power(t1,3) - 6*Power(t1,2)*t2 + 2*t1*Power(t2,2) - xihat))/ (6*Power(t1,3) - 6*Power(t1,2)*t2 + 2*t1*Power(t2,2) - 4*xihat) + 6*Power(a1,2)*(-3*t1 + t2)*(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat) + (3*xihat*(9*Power(t1,6) - 18*Power(t1,5)*t2 + 15*Power(t1,4)*Power(t2,2) - 6*Power(t1,3)*(Power(t2,3) - 7*xihat) + 14*t1*Power(t2,2)*xihat + 4*Power(xihat,2) + Power(t1,2)*(Power(t2,4) - 42*t2*xihat)))/ (2.*(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) - 2*xihat)* Power(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat,2)) + (a1*Power(-3*t1 + t2,2)*(18*Power(t1,6) - 36*Power(t1,5)*t2 + 30*Power(t1,4)*Power(t2,2) + t1*Power(t2,2)*xihat + 8*Power(xihat,2) + 3*Power(t1,3)*(-4*Power(t2,3) + xihat) + Power(t1,2)*(2*Power(t2,4) - 3*t2*xihat)))/ ((3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) - 2*xihat)* (3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat))))/Power(E,a1*Power(-3*t1 + t2,2))\ + (36*A1*W0*((xihat*(Power(t2,6) - 2*Power(t2,3)*Power(-3*t1 + t2,3) + Power(-3*t1 + t2,6) + 126*Power(t2,3)*xihat - 126*Power(-3*t1 + t2,3)*xihat + 324*Power(xihat,2)))/108. + (a1*Power(-3*t1 + t2,2)*(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat)* (18*Power(t1,6) - 36*Power(t1,5)*t2 + 30*Power(t1,4)*Power(t2,2) + t1*Power(t2,2)*xihat + 8*Power(xihat,2) + 3*Power(t1,3)*(-4*Power(t2,3) + xihat) + Power(t1,2)*(2*Power(t2,4) - 3*t2*xihat)))/4.)*Cos(a1*r1))/ (Power(E,(a1*Power(-3*t1 + t2,2))/2.)*(Power(3*t1 - t2,3) + Power(t2,3) - 18*xihat)* Power(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat,2)) + (36*A2*Power(E,(a2*t1*(3*t1 - 2*t2))/2.)*W0* ((xihat*(Power(t2,6) - 2*Power(t2,3)*Power(-3*t1 + t2,3) + Power(-3*t1 + t2,6) + 126*Power(t2,3)*xihat - 126*Power(-3*t1 + t2,3)*xihat + 324*Power(xihat,2)))/108. - (a2*t1*(3*t1 - 2*t2)*(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat)* (18*Power(t1,6) - 36*Power(t1,5)*t2 + 30*Power(t1,4)*Power(t2,2) + t1*Power(t2,2)*xihat + 8*Power(xihat,2) + 3*Power(t1,3)*(-4*Power(t2,3) + xihat) + Power(t1,2)*(2*Power(t2,4) - 3*t2*xihat)))/4.)*Cos(a2*r2))/ ((Power(3*t1 - t2,3) + Power(t2,3) - 18*xihat)* Power(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat,2)) + A1*A2*Power(E,(a2*t1*(3*t1 - 2*t2) - a1*Power(-3*t1 + t2,2))/2.)* (-((a1*a2*t1*(3*t1 - 2*t2)*Power(-3*t1 + t2,2)* (6*Power(t1,3) - 6*Power(t1,2)*t2 + 2*t1*Power(t2,2) - xihat))/ (3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) - 2*xihat)) - 4*a1*a2*(-3*t1 + t2)*(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat) + (3*xihat*(9*Power(t1,6) - 18*Power(t1,5)*t2 + 15*Power(t1,4)*Power(t2,2) - 6*Power(t1,3)*(Power(t2,3) - 7*xihat) + 14*t1*Power(t2,2)*xihat + 4*Power(xihat,2) + Power(t1,2)*(Power(t2,4) - 42*t2*xihat)))/ ((3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) - 2*xihat)* Power(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat,2)) + ((a1*Power(-3*t1 + t2,2) + a2*t1*(-3*t1 + 2*t2))* (18*Power(t1,6) - 36*Power(t1,5)*t2 + 30*Power(t1,4)*Power(t2,2) + t1*Power(t2,2)*xihat + 8*Power(xihat,2) + 3*Power(t1,3)*(-4*Power(t2,3) + xihat) + Power(t1,2)*(2*Power(t2,4) - 3*t2*xihat)))/ ((3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) - 2*xihat)* (3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat)))*Cos(a1*r1 - a2*r2)))/ Power(3*Power(t1,3) - 3*Power(t1,2)*t2 + t1*Power(t2,2) + xihat,2) ) V = Vc.real return V ############################################################################################# ############################################################################################# ############################# variety of functions defined above ############################ ########################### negf is the one that is returned below ########################## ############################################################################################# ############################################################################################# ########################### parameters selected below ####################################### ############################################################################################# pi=3.1415926535897932385 ################################ parameters _racetrack ###################################### Kcs = -3.50656 k111 = 1.0 A1 = 1 A2 = 1.1 a1 = pi/11 a2 = pi/12 W0 = -0.812 gstring = 0.14306151645207438 chiiP11169 = -186.0 params_racetrack = np.array([Kcs,k111,chiiP11169,W0,A1,a1,A2,a2,gstring]) content_racetrack= '\n' + r'''The parameters used in this example are \begin{align} g_{string} & = '''+ str(gstring) + r''' ~;~~ \chi = '''+ str(chiiP11169) + r'''~;~~ W_0 = ''' + str(W0) + r'''\nonumber \\ A_1 &= ''' + str(A1) + r'''~;~~ a_1 = ''' + str(a1) + r'''~;~~ A_2 = ''' + str(A2) + r'''\nonumber \\ a_2 &= ''' + str(a2) + r'''~;~~ K_{cs} = ''' + str(Kcs) + r'''~;~~ k_{111} = ''' + str(k111) + r'''~;~~ \end{align} ''' ############################### parameters wobbly gaussian ################################## params_wobbly_gaussian=np.array([]) content_wobbly_gaussian= '\n' + r'''The Wobbly Gaussian''' ############################################################################################# ############################################################################################# ############################################################################################# ##### select function below (by alterining the _extension) -- negf is the POTENTIAL ######### ##### params is the parameters if there are ##### content2 is the latex content to accompany the plot ##### identify the variables in the main programme ########################################## params=params_wobbly_gaussian content2=content_wobbly_gaussian names=names_racetrack # choose your names here def negf(xx): # choose the fitness f, aka negative of potential - negf is negative fitness return negf_wobbly_guassian(xx) #return negf_racetrack(params,xx) def potential(xx): # potential return negf_wobbly_guassian(xx) #return negf_racetrack(params,xx) ############################################################################################# #################### don't edit below unless you know what you are doing #################### ############################################################################################# def f(xx): return -negf(xx) def func(i,j,qq,xx,yy): # dummy function for plotting that calls the real function with x,y inserted in the right entries. xx are the coords -- this doesn't work for me but I leave it here anyway np.put(qq,[i,j],[xx,yy]) return f(qq) def isNaN(num): return num != num def close(aa,bb): top=aa-bb bot=aa+bb print (top,bot) vec=top for i in xrange(0,len(top)): print (abs(top[i]/bot[i])) if abs(top[i]/bot[i]) > 0.001: return False return True