Why this function always results in a wave with value 0?

Share icon

Thu, 11/09/2023 - 08:30 pm

make/d/o/n=20 twave=-10^(6)+2*10^(6)/19*p
make/d/o/n=40 zwave=10^(9)*p
make/d/o/n=(20,40)/c yw

function/c f1(inp)
variable inp
nvar ini,a,b,c,m,hbar
nvar t,z
variable inp0=sqrt(2*m*0.1),v0=sqrt(2*0.1/m)
variable/c temp1=cmplx(0,-(inp)^2/2/m*t/hbar)
variable/c temp2=cmplx(0,(inp)*z/hbar)
variable phi0=pi/10
variable/c temp3=cmplx(0,ini*phi0)
variable temp4=abs(2*a)
return besselj(ini,temp4)*exp(-(inp-inp0-2.40747*10^(-6)/v0*ini)^2/4/b^2)/sqrt(sqrt(2*pi)/b)*exp(temp1)*exp(temp2)*exp(temp3)
end

function/c f10(tt,zz)
variable tt,zz
nvar ini,a,b,c,m,hbar
variable/g t,z
t=tt
z=zz
variable/c cir=integrate1d(f1,-100,100,0)
return cir
end


function ffinal(twave,zwave,yw)
wave twave,zwave
wave/c yw
variable/g ini,a=2,c=1,m=0.511,hbar=1
variable/g b=m*10^(-6)/sqrt(2*m*0.1)
variable i,j,ii
yw=0
for(i=0;i<20;i++)
  for(j=0;j<40;j++)
    for(ii=-5;ii<6;ii++)
      ini=ii
      yw[i][j]+=f10(twave[i],zwave[j])
    endfor
  endfor
endfor
end

I want to simulate PINEM electron density probability as a function of time and position after interaction with near field. But the function ffinal(twave,zwave,yw) always results in a zero yw wave. Basically it is the time-dependent momentum space electron wavefunction Fourier transformed to position space.What I want to get is shown in the attachment. It should be magsqr(yw). But yw in my code always equals zero. What's wrong with my code? Many thanks!

 PINEM electron density probability as a function of time and position after interaction with near field

I'm betting you copied source from another program (perhaps MatLab)?

In any event, it is helpful to both use the debugger and to rewrite complicated expressions in pieces so you can see what the numeric results of each piece is.

In my rewritten f1 code, argReal = -984739315945430, and exp( -984739315945430) is identically 0. It isn't going to matter what you multiply it with, the result will be 0.

I think you have one or more typos in f1()

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3             // Use modern global access method and strict wave access
#pragma DefaultTab={3,20,4}     // Set default tab width in Igor Pro 9 and later

Macro TryIt()
    make/d/o/n=20 twave=-10^(6)+2*10^(6)/19*p
    make/d/o/n=40 zwave=10^(9)*p
    make/d/o/n=(20,40)/c yw
    ffinal(twave,zwave,yw)
End

function/c f1(inp)
    variable inp
    nvar ini,a,b,c,m,hbar
    nvar t,z
    variable inp0=sqrt(2*m*0.1),v0=sqrt(2*0.1/m)
    variable/c temp1=cmplx(0,-(inp)^2/2/m*t/hbar)
    variable/c temp2=cmplx(0,(inp)*z/hbar)
    variable phi0=pi/10
    variable/c temp3=cmplx(0,ini*phi0)
    variable temp4=abs(2*a)
    Variable aReal=besselj(ini,temp4)
   
    Variable/c expt1 = exp(temp1)
    Variable/c expt2 = exp(temp2)
    Variable/C expt3 = exp(temp3)
    Variable/C prod = expt1*expt2*expt3
    Variable argReal = -(inp-inp0-2.40747*10^-6/v0*ini)^2/4/b^2
    Variable denomReal = sqrt(sqrt(2*pi)/b)
    Variable/C result = aReal*exp(argReal)/denomReal*prod
    return result
End


function/c f10(tt,zz)
    variable tt,zz
    nvar ini,a,b,c,m,hbar
    variable/g t,z
    t=tt
    z=zz
    variable/c cir=integrate1d(f1,-100,100,0)
    return cir
end


function ffinal(twave,zwave,yw)
    wave twave,zwave
    wave/c yw
    variable/g ini,a=2,c=1,m=0.511,hbar=1
    variable/g b=m*10^(-6)/sqrt(2*m*0.1)
    variable i,j,ii
    yw=0
    for(i=0;i<20;i++)
        for(j=0;j<40;j++)
            for(ii=-5;ii<6;ii++)
                ini=ii
                Variable/C addThis = f10(twave[i],zwave[j])
                yw[i][j]+=addThis
            endfor
        endfor
    endfor
end

November 9, 2023 at 09:53 pm - Permalink

Hello Yuan,

It might help if you posted an image of the original equation that you are trying to evaluate together with the values of the relevant parameters.  As Jim indicated above, the numerical factors that you are evaluating do not lead to a useful result.  This is similar conceptually to implementing the sinc(x) function by directly computing sin(x)/x at x=0.

 

A.G.

November 10, 2023 at 12:20 pm - Permalink