From e72af82c3ed16b81cdb5043c58abbdbb3cf02102 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Mon, 4 Oct 2021 11:12:34 -0400 Subject: Initial commit --- src/anyeq.jl | 949 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 949 insertions(+) create mode 100644 src/anyeq.jl (limited to 'src/anyeq.jl') diff --git a/src/anyeq.jl b/src/anyeq.jl new file mode 100644 index 0000000..8459cb0 --- /dev/null +++ b/src/anyeq.jl @@ -0,0 +1,949 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# interpolation +@everywhere mutable struct Anyeq_approx + aK::Float64 + bK::Float64 + gK::Float64 + aL1::Float64 + bL1::Float64 + aL2::Float64 + bL2::Float64 + gL2::Float64 + aL3::Float64 + bL3::Float64 + gL3::Float64 +end + +# compute energy for a given rho +# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) +function anyeq_energy(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) + # initialize vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # compute initial guess from medeq + rhos=Array{Float64}(undef,nlrho) + for j in 0:nlrho-1 + rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) + end + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + u0=u0s[nlrho] + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + + @printf("% .15e % .15e\n",E,error) +end + +# compute energy as a function of rho +function anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) + # initialize vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # compute initial guess from medeq + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + + # save result from each task + es=Array{Float64,1}(undef,length(rhos)) + err=Array{Float64,1}(undef,length(rhos)) + + ## spawn workers + # number of workers + nw=nworkers() + # split jobs among workers + work=Array{Array{Int64,1},1}(undef,nw) + # init empty arrays + for p in 1:nw + work[p]=zeros(0) + end + for j in 1:length(rhos) + append!(work[(j-1)%nw+1],j) + end + + count=0 + # for each worker + @sync for p in 1:nw + # for each task + @async for j in work[p] + count=count+1 + if count>=nw + progress(count,length(rhos),10000) + end + # run the task + (u,es[j],err[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + end + end + + for j in 1:length(rhos) + @printf("% .15e % .15e % .15e\n",rhos[j],es[j],err[j]) + end +end + +# compute energy as a function of rho +# initialize with previous rho +function anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) + # initialize vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # compute initial guess from medeq + u0s=anyeq_init_medeq([rhos[1]],N,J,k,a0,v,maxiter,tolerance) + u=u0s[1] + + for j in 1:length(rhos) + progress(j,length(rhos),10000) + # run the task + # init Newton from previous rho + (u,E,error)=anyeq_hatu(u,P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + @printf("% .15e % .15e % .15e\n",rhos[j],E,error) + # abort when the error gets too big + if error>tolerance + break + end + end +end +# compute energy as a function of rho +# initialize with next rho +function anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) + # initialize vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # compute initial guess from medeq + u0s=anyeq_init_medeq([rhos[length(rhos)]],N,J,k,a0,v,maxiter,tolerance) + u=u0s[1] + + for j in 1:length(rhos) + progress(j,length(rhos),10000) + # run the task + # init Newton from previous rho + (u,E,error)=anyeq_hatu(u,P,N,J,rhos[length(rhos)+1-j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + @printf("% .15e % .15e % .15e\n",rhos[length(rhos)+1-j],real(E),error) + # abort when the error gets too big + if error>tolerance + break + end + end +end + +# compute u(k) +# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) +function anyeq_uk(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,approx,savefile) + # init vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + # compute initial guess from medeq + rhos=Array{Float64}(undef,nlrho) + for j in 0:nlrho-1 + rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) + end + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + u0=u0s[nlrho] + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + + for zeta in 0:J-1 + for j in 1:N + # order k's in increasing order + @printf("% .15e % .15e\n",k[(J-1-zeta)*N+j],u[(J-1-zeta)*N+j]) + end + end +end + +# compute u(x) +# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) +function anyeq_ux(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile) + # init vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # compute initial guess from medeq + rhos=Array{Float64}(undef,nlrho) + for j in 0:nlrho-1 + rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) + end + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + u0=u0s[nlrho] + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + + for i in 1:nx + x=xmin+(xmax-xmin)*i/nx + ux=0. + for zeta in 0:J-1 + for j in 1:N + ux+=(taus[zeta+2]-taus[zeta+1])/(16*pi*x)*weights[2][j]*cos(pi*weights[1][j]/2)*(1+k[zeta*N+j])^2*k[zeta*N+j]*u[zeta*N+j]*sin(k[zeta*N+j]*x) + end + end + @printf("% .15e % .15e % .15e\n",x,real(ux),imag(ux)) + end +end + +# compute condensate fraction for a given rho +# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) +function anyeq_condensate_fraction(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) + # initialize vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # compute initial guess from medeq + rhos=Array{Float64}(undef,nlrho) + for j in 0:nlrho-1 + rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) + end + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + u0=u0s[nlrho] + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + + # compute eta + eta=anyeq_eta(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + + @printf("% .15e % .15e\n",eta,error) +end + +# condensate fraction as a function of rho +function anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) + ## spawn workers + # number of workers + nw=nworkers() + # split jobs among workers + work=Array{Array{Int64,1},1}(undef,nw) + # init empty arrays + for p in 1:nw + work[p]=zeros(0) + end + for j in 1:length(rhos) + append!(work[(j-1)%nw+1],j) + end + + # initialize vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # compute initial guess from medeq + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + + # compute u + us=Array{Array{Float64,1}}(undef,length(rhos)) + errs=Array{Float64,1}(undef,length(rhos)) + count=0 + # for each worker + @sync for p in 1:nw + # for each task + @async for j in work[p] + count=count+1 + if count>=nw + progress(count,length(rhos),10000) + end + # run the task + (us[j],E,errs[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + end + end + + # compute eta + etas=Array{Float64}(undef,length(rhos)) + count=0 + # for each worker + @sync for p in 1:nw + # for each task + @async for j in work[p] + count=count+1 + if count>=nw + progress(count,length(rhos),10000) + end + # run the task + etas[j]=anyeq_eta(us[j],P,N,J,rhos[j],weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + end + end + + for j in 1:length(rhos) + @printf("% .15e % .15e % .15e\n",rhos[j],etas[j],errs[j]) + end +end + +# compute the momentum distribution for a given rho +# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) +function anyeq_momentum_distribution(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) + # initialize vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # compute initial guess from medeq + rhos=Array{Float64}(undef,nlrho) + for j in 0:nlrho-1 + rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) + end + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + u0=u0s[nlrho] + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + + # compute M + M=anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + + for zeta in 0:J-1 + for j in 1:N + # order k's in increasing order + @printf("% .15e % .15e\n",k[(J-1-zeta)*N+j],M[(J-1-zeta)*N+j]) + end + end +end + +# 2 point correlation function +function anyeq_2pt_correlation(minlrho,nlrho,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile) + # init vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # compute initial guess from medeq + rhos=Array{Float64}(undef,nlrho) + for j in 0:nlrho-1 + rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) + end + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + u0=u0s[nlrho] + + # compute u and some useful integrals + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + for i in 1:nx + x=xmin+(xmax-xmin)*i/nx + C2=anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK,G) + @printf("% .15e % .15e\n",x,C2) + end +end + +# compute Abar +function anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + if approx.bL3==0. + return [] + end + + out=Array{Array{Float64,5}}(undef,J*N) + + ## spawn workers + # number of workers + nw=nworkers() + # split jobs among workers + work=Array{Array{Int64,1},1}(undef,nw) + # init empty arrays + for p in 1:nw + work[p]=zeros(0) + end + for j in 1:N*J + append!(work[(j-1)%nw+1],j) + end + + count=0 + # for each worker + @sync for p in 1:nw + # for each task + @async for j in work[p] + count=count+1 + if count>=nw + progress(count,N*J,10000) + end + # run the task + out[j]=remotecall_fetch(barAmat,workers()[p],k[j],weights,taus,T,P,N,J,2,2,2,2,2) + end + end + + return out +end + +# initialize computation +@everywhere function anyeq_init(taus,P,N,J,v) + # Gauss-Legendre weights + weights=gausslegendre(N) + + # initialize vectors V,k + V=Array{Float64}(undef,J*N) + k=Array{Float64}(undef,J*N) + for zeta in 0:J-1 + for j in 1:N + xj=weights[1][j] + # set kj + k[zeta*N+j]=(2+(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)-(taus[zeta+2]+taus[zeta+1]))/(2-(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)+taus[zeta+2]+taus[zeta+1]) + # set v + V[zeta*N+j]=v(k[zeta*N+j]) + end + end + # potential at 0 + V0=v(0) + + # initialize matrix A + T=chebyshev_polynomials(P) + A=Amat(k,weights,taus,T,P,N,J,2,2) + + # compute Upsilon + # Upsilonmat does not use splines, so increase precision + weights_plus=gausslegendre(N*J) + Upsilon=Upsilonmat(k,v,weights_plus) + Upsilon0=Upsilon0mat(k,v,weights_plus) + + return(weights,T,k,V,V0,A,Upsilon,Upsilon0) +end + +# compute initial guess from medeq +@everywhere function anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + us_medeq=Array{Array{Float64,1}}(undef,length(rhos)) + u0s=Array{Array{Float64,1}}(undef,length(rhos)) + + weights_medeq=gausslegendre(N*J) + + (us_medeq[1],E,err)=easyeq_hatu(easyeq_init_u(a0,J*N,weights_medeq),J*N,rhos[1],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.)) + u0s[1]=easyeq_to_anyeq(us_medeq[1],weights_medeq,k,N,J) + if err>tolerance + print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[1],"\n") + end + + for j in 2:length(rhos) + (us_medeq[j],E,err)=easyeq_hatu(us_medeq[j-1],J*N,rhos[j],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.)) + u0s[j]=easyeq_to_anyeq(us_medeq[j],weights_medeq,k,N,J) + + if err>tolerance + print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[j],"\n") + end + end + + return u0s +end + +# interpolate the solution of medeq to an input for anyeq +@everywhere function easyeq_to_anyeq(u_simple,weights,k,N,J) + # reorder u_simple, which is evaluated at (1-x_j)/(1+x_j) with x_j\in[-1,1] + u_s=zeros(Float64,length(u_simple)) + k_s=Array{Float64}(undef,length(u_simple)) + for j in 1:length(u_simple) + xj=weights[1][j] + k_s[length(u_simple)-j+1]=(1-xj)/(1+xj) + u_s[length(u_simple)-j+1]=u_simple[j] + end + + # initialize U + u=zeros(Float64,J*N) + for zeta in 0:J-1 + for j in 1:N + u[zeta*N+j]=linear_interpolation(k[zeta*N+j],k_s,u_s) + end + end + + return u +end + + +# compute u using chebyshev expansions +@everywhere function anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + # init + # rescale by rho (that's how u is defined) + u=rho*u0 + + # quantify relative error + error=-1. + + # run Newton algorithm + for i in 1:maxiter-1 + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + new=u-inv(anyeq_DXi(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_Xi(u,X,Y) + + error=norm(new-u)/norm(u) + if(error0 + dV[i]=sin(k[i]*x)/(k[i]*x)*hann(k[i],windowL) + else + dV[i]=hann(k[i],windowL) + end + end + dV0=1. + + # compute dUpsilon + # Upsilonmat does not use splines, so increase precision + weights_plus=gausslegendre(N*J) + dUpsilon=Upsilonmat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus) + dUpsilon0=Upsilon0mat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus) + + du=-inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_dXidv(x,rho*u,rho,k,taus,dV,dV0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) + # rescale rho + du=du/rho + + C2=rho^2*(1-integrate_f_chebyshev(s->1.,u.*dV+V.*du,k,taus,weights,N,J)) + + return C2 +end + +# derivative of Xi with respect to v in the direction sin(kx)/kx +@everywhere function anyeq_dXidv(x,U,rho,k,taus,dv,dv0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) + # Chebyshev expansion of U + FU=chebyshev(U,taus,weights,P,N,J,2) + + dS=dv-conv_v_chebyshev(U,dUpsilon,k,taus,weights,N,J)/rho + dE=dv0-avg_v_chebyshev(U,dUpsilon0,k,taus,weights,N,J)/rho + + UU=conv_chebyshev(FU,FU,A) + + dII=zeros(Float64,N*J) + if approx.gL2!=0. + if approx.bL2!=0. + dII+=approx.gL2*approx.bL2*conv_chebyshev(FU,chebyshev(U.*dS,taus,weights,P,N,J,2),A)/rho + end + if approx.bL2!=1. + dII+=approx.gL2*(1-approx,bL2)*dE/rho*UU + end + end + + dJJ=zeros(Float64,J*N) + if approx.gL3!=0. + if approx.bL3!=0. + dFS=chebyshev(dS,taus,weights,P,N,J,2) + dJJ+=approx.gL3*approx.bL3*double_conv_S_chebyshev(FU,FU,FU,FU,dFS,Abar) + end + if approx.bL3!=1. + dJJ=approx.gL3*(1-approx.bL3)*dE*(UU/rho).^2 + end + end + + dG=zeros(Float64,N*J) + if approx.aK!=0. && approx.gK!=0. + if approx.bK!=0. + dG+=approx.aK*approx.gK*approx.bK*conv_chebyshev(FU,chebyshev(2*dS.*U,taus,weights,P,N,J,2),A)/rho + end + if approx.bK!=1. + dG+=approx.aK*approx.gK*(1-approx.bK)*2*dE*UU/rho + end + end + if approx.aL1!=0. + if approx.bL1!=0. + dG-=approx.aL1*approx.bL1*conv_chebyshev(FU,chebyshev(dS.*(U.^2),taus,weights,P,N,J,2),A)/rho + end + if approx.bL1!=1. + dG-=approx.aL1*(1-approx.bL1)*dE/rho*conv_chebyshev(FU,chebyshev((U.^2),taus,weights,P,N,J,2),A) + end + end + if approx.aL2!=0. && approx.gL2!=0. + dG+=approx.aL2*conv_chebyshev(FU,chebyshev(2*dII.*U,taus,weights,P,N,J,2),A)/rho + end + if approx.aL3!=0. && approx.gL3!=0. + dG-=approx.aL3*conv_chebyshev(FU,chebyshev(dJJ/2,taus,weights,P,N,J,2),A)/rho + end + + dsK=zeros(Float64,N*J) + if approx.gK!=0. + if approx.bK!=0. + dsK+=approx.gK*approx.bK*dS + end + if approx.bK!=1. + dsK+=approx.gK*(1-approx.bK)*dE*ones(N*J) + end + end + dsL1=zeros(Float64,N*J) + if approx.bL1!=0. + dsL1+=approx.bL1*dS + end + if approx.bL1!=1. + dsL1+=(1-approx.bL1)*dE*ones(N*J) + end + + dX=(dsK-dsL1+dII)./sL1-X./sL1.*dsL1 + dY=(dS-dsL1+dG+dJJ/2)./sL1-Y./sL1.*dsL1 + + out=((Y.+1).*dX./(X.+1)-dY)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1)./(2*(X.+1).^3).*(2*(Y.+1)./(X.+1).*dX-dY).*dotdPhi((Y.+1)./(X.+1).^2) + return out +end -- cgit v1.2.3-54-g00ecf