![](/profiles/wavemetrics/themes/wavemetrics/logo.png)
Plot COMSOL multiphysics PINEM simulation
![](/sites/default/files/styles/thumbnail/public/Yuan%20Fang/profile-images/2020-11/01-1589457975-834716.jpg?itok=rYibQIgu)
Yuan Fang
I wrote a code to plot photon-induced near field electron microscopy (PINEM) coupling constant simulation using COMSOL multiphysics. The code is attached below:
function readcomsol(path1,fname1)
string path1,fname1
variable a,i,t,j,k
string str,fword,nline,nbnds
Open/R a path1 + fname1
for(i=0;i<4;i++)
FReadLine a, str
endfor
FReadLine a, str
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==3)
nline=fword
endif
endif
i=i+1
while(t<3)
make/d/n=(str2num(nline),4) w
for(i=0;i<3;i++)
FReadLine a, str
endfor
for(i=0;i<str2num(nline);i++)
FReadLine a, str
k=0
t=0
do
fword=stringfromlist(k,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
w[i][0]=str2num(fword)
endif
if(t==2)
w[i][1]=str2num(fword)
endif
if(t==3)
w[i][2]=str2num(fword)
endif
if(t==4)
w[i][3]=str2num(fword)
endif
endif
k=k+1
while(t<4)
endfor
close a
duplicate/o/rmd=[0,str2num(nline)-1][0,0] w,twave
make/n=(str2num(nline))/d/o xwave=twave
duplicate/o/rmd=[0,str2num(nline)-1][1,1] w,twave
make/n=(str2num(nline))/d/o ywave=twave
duplicate/o/rmd=[0,str2num(nline)-1][2,2] w,twave
make/n=(str2num(nline))/d/o zwave=twave
duplicate/o/rmd=[0,str2num(nline)-1][3,3] w,twave
make/n=(str2num(nline))/d/o fwave=twave
end
Function plotcomsol(pixels,w,plane)
variable pixels
wave w
string plane
Variable xMin, xMax, yMin, yMax
Make /O /N=(DimSize(w,0)) tempWave=0
if(stringmatch(plane,"xy"))
tempWave = w[p][0]
WaveStats /Q tempWave
xMin = V_min
xMax = V_max
tempWave = w[p][1]
WaveStats /Q tempWave
yMin = V_min
yMax = V_max
endif
if(stringmatch(plane,"yz"))
tempWave = w[p][1]
WaveStats /Q tempWave
xMin = V_min
xMax = V_max
tempWave = w[p][2]
WaveStats /Q tempWave
yMin = V_min
yMax = V_max
endif
if(stringmatch(plane,"xz"))
tempWave = w[p][0]
WaveStats /Q tempWave
xMin = V_min
xMax = V_max
tempWave = w[p][2]
WaveStats /Q tempWave
yMin = V_min
yMax = V_max
endif
// Determine the final size for the map
Variable ratio = (xMax - xMin)/(yMax - yMin)
Variable mapSizeX, mapSizeY, delta,deltai
if (ratio > 1)
mapSizeX = pixels
mapSizeY = floor(pixels/ratio)
delta = (xMax - xMin)/(pixels - 1)
else
mapSizeX = floor(pixels*ratio)
mapSizeY = pixels
delta = (yMax - yMin)/(pixels - 1)
endif
make/d/o /N=(mapSizeX, mapSizeY) map
SetScale /P x xMin, delta, map
SetScale /P y yMin, delta, map
map=0
duplicate/o map,BlueJacobian
Variable i,j,pee, qew,u,v
for (j=0; j<DimSize(w,0); j+=1)
if(stringmatch(plane,"xy"))
pee = round((w[j][0] - DimOffset(map,0))/DimDelta(map,0))
qew = round((w[j][1] - DimOffset(map,1))/DimDelta(map,1))
endif
if(stringmatch(plane,"yz"))
pee = round((w[j][1] - DimOffset(map,0))/DimDelta(map,0))
qew = round((w[j][2] - DimOffset(map,1))/DimDelta(map,1))
endif
if(stringmatch(plane,"xz"))
pee = round((w[j][0] - DimOffset(map,0))/DimDelta(map,0))
qew = round((w[j][2] - DimOffset(map,1))/DimDelta(map,1))
endif
if (pee >= 0 && pee < DimSize(map,0) && qew >= 0 && qew < DimSize(map,1))
map[pee][qew] += w[j][3]
BlueJacobian[pee][qew] += 1
endif
endfor
for(u=0;u<DimSize(map,0);u+=1)
for(v=0;v<DimSize(map,1);v+=1)
map[u][v] /= max(1,BlueJacobian[u][v])
endfor
endfor
display
appendimage map
ModifyImage map ctab= {*,*,ColdWarm,0}
for(i=0;i<dimsize(map,0);i++)
for(j=0;j<dimsize(map,1);j++)
if(i==0&&j==0)
map[i][j]=map[i][j]==0 ? (map[i+1][j]+map[i][j+1])/2 : map[i][j]
elseif(i==0&&j==(dimsize(map,1)-1))
map[i][j]=map[i][j]==0 ? (map[i][j-1]+map[i+1][j])/2 : map[i][j]
elseif(i==(dimsize(map,0)-1)&&j==0)
map[i][j]=map[i][j]==0 ? (map[i][j+1]+map[i-1][j])/2 : map[i][j]
elseif(i==(dimsize(map,0)-1)&&j==(dimsize(map,1)-1))
map[i][j]=map[i][j]==0 ? (map[i-1][j]+map[i][j-1])/2 : map[i][j]
elseif(i==0&&j!=(dimsize(map,1)-1)&&j!=0)
map[i][j]=map[i][j]==0 ? (map[i][j-1]+map[i+1][j]+map[i][j+1])/3 : map[i][j]
elseif(j==0&&i!=(dimsize(map,1)-1)&&i!=0)
map[i][j]=map[i][j]==0 ? (map[i-1][j]+map[i+1][j]+map[i][j+1])/3 : map[i][j]
elseif(i==(dimsize(map,0)-1)&&j!=(dimsize(map,1)-1)&&j!=0)
map[i][j]=map[i][j]==0 ? (map[i][j-1]+map[i-1][j]+map[i][j+1])/3 : map[i][j]
elseif(j==(dimsize(map,0)-1)&&i!=(dimsize(map,1)-1)&&i!=0)
map[i][j]=map[i][j]==0 ? (map[i][j-1]+map[i-1][j]+map[i+1][j])/3 : map[i][j]
//else
// map[i][j]=map[i][j]==0 ? (map[i-1][j]+map[i][j-1]+map[i+1][j]+map[i][j+1])/4 : map[i][j]
endif
endfor
endfor
// make/d/o/n=(mapsizex/2,mapsizey/2) mapc
// SetScale /P x xMin, 2*delta, mapc
// SetScale /P y yMin, 2*delta, mapc
// mapc=0
// for(i=0;i<dimsize(mapc,0);i++)
// for(j=0;j<dimsize(mapc,1);j++)
// mapc[i][j]+=map[2*i][2*j]
// mapc[i][j]+=map[2*i+1][2*j]
// mapc[i][j]+=map[2*i][2*j+1]
// mapc[i][j]+=map[2*i+1][2*j+1]
// endfor
// endfor
// display
// appendimage mapc
// ModifyImage mapc ctab= {*,*,ColdWarm,0}
End
string path1,fname1
variable a,i,t,j,k
string str,fword,nline,nbnds
Open/R a path1 + fname1
for(i=0;i<4;i++)
FReadLine a, str
endfor
FReadLine a, str
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==3)
nline=fword
endif
endif
i=i+1
while(t<3)
make/d/n=(str2num(nline),4) w
for(i=0;i<3;i++)
FReadLine a, str
endfor
for(i=0;i<str2num(nline);i++)
FReadLine a, str
k=0
t=0
do
fword=stringfromlist(k,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
w[i][0]=str2num(fword)
endif
if(t==2)
w[i][1]=str2num(fword)
endif
if(t==3)
w[i][2]=str2num(fword)
endif
if(t==4)
w[i][3]=str2num(fword)
endif
endif
k=k+1
while(t<4)
endfor
close a
duplicate/o/rmd=[0,str2num(nline)-1][0,0] w,twave
make/n=(str2num(nline))/d/o xwave=twave
duplicate/o/rmd=[0,str2num(nline)-1][1,1] w,twave
make/n=(str2num(nline))/d/o ywave=twave
duplicate/o/rmd=[0,str2num(nline)-1][2,2] w,twave
make/n=(str2num(nline))/d/o zwave=twave
duplicate/o/rmd=[0,str2num(nline)-1][3,3] w,twave
make/n=(str2num(nline))/d/o fwave=twave
end
Function plotcomsol(pixels,w,plane)
variable pixels
wave w
string plane
Variable xMin, xMax, yMin, yMax
Make /O /N=(DimSize(w,0)) tempWave=0
if(stringmatch(plane,"xy"))
tempWave = w[p][0]
WaveStats /Q tempWave
xMin = V_min
xMax = V_max
tempWave = w[p][1]
WaveStats /Q tempWave
yMin = V_min
yMax = V_max
endif
if(stringmatch(plane,"yz"))
tempWave = w[p][1]
WaveStats /Q tempWave
xMin = V_min
xMax = V_max
tempWave = w[p][2]
WaveStats /Q tempWave
yMin = V_min
yMax = V_max
endif
if(stringmatch(plane,"xz"))
tempWave = w[p][0]
WaveStats /Q tempWave
xMin = V_min
xMax = V_max
tempWave = w[p][2]
WaveStats /Q tempWave
yMin = V_min
yMax = V_max
endif
// Determine the final size for the map
Variable ratio = (xMax - xMin)/(yMax - yMin)
Variable mapSizeX, mapSizeY, delta,deltai
if (ratio > 1)
mapSizeX = pixels
mapSizeY = floor(pixels/ratio)
delta = (xMax - xMin)/(pixels - 1)
else
mapSizeX = floor(pixels*ratio)
mapSizeY = pixels
delta = (yMax - yMin)/(pixels - 1)
endif
make/d/o /N=(mapSizeX, mapSizeY) map
SetScale /P x xMin, delta, map
SetScale /P y yMin, delta, map
map=0
duplicate/o map,BlueJacobian
Variable i,j,pee, qew,u,v
for (j=0; j<DimSize(w,0); j+=1)
if(stringmatch(plane,"xy"))
pee = round((w[j][0] - DimOffset(map,0))/DimDelta(map,0))
qew = round((w[j][1] - DimOffset(map,1))/DimDelta(map,1))
endif
if(stringmatch(plane,"yz"))
pee = round((w[j][1] - DimOffset(map,0))/DimDelta(map,0))
qew = round((w[j][2] - DimOffset(map,1))/DimDelta(map,1))
endif
if(stringmatch(plane,"xz"))
pee = round((w[j][0] - DimOffset(map,0))/DimDelta(map,0))
qew = round((w[j][2] - DimOffset(map,1))/DimDelta(map,1))
endif
if (pee >= 0 && pee < DimSize(map,0) && qew >= 0 && qew < DimSize(map,1))
map[pee][qew] += w[j][3]
BlueJacobian[pee][qew] += 1
endif
endfor
for(u=0;u<DimSize(map,0);u+=1)
for(v=0;v<DimSize(map,1);v+=1)
map[u][v] /= max(1,BlueJacobian[u][v])
endfor
endfor
display
appendimage map
ModifyImage map ctab= {*,*,ColdWarm,0}
for(i=0;i<dimsize(map,0);i++)
for(j=0;j<dimsize(map,1);j++)
if(i==0&&j==0)
map[i][j]=map[i][j]==0 ? (map[i+1][j]+map[i][j+1])/2 : map[i][j]
elseif(i==0&&j==(dimsize(map,1)-1))
map[i][j]=map[i][j]==0 ? (map[i][j-1]+map[i+1][j])/2 : map[i][j]
elseif(i==(dimsize(map,0)-1)&&j==0)
map[i][j]=map[i][j]==0 ? (map[i][j+1]+map[i-1][j])/2 : map[i][j]
elseif(i==(dimsize(map,0)-1)&&j==(dimsize(map,1)-1))
map[i][j]=map[i][j]==0 ? (map[i-1][j]+map[i][j-1])/2 : map[i][j]
elseif(i==0&&j!=(dimsize(map,1)-1)&&j!=0)
map[i][j]=map[i][j]==0 ? (map[i][j-1]+map[i+1][j]+map[i][j+1])/3 : map[i][j]
elseif(j==0&&i!=(dimsize(map,1)-1)&&i!=0)
map[i][j]=map[i][j]==0 ? (map[i-1][j]+map[i+1][j]+map[i][j+1])/3 : map[i][j]
elseif(i==(dimsize(map,0)-1)&&j!=(dimsize(map,1)-1)&&j!=0)
map[i][j]=map[i][j]==0 ? (map[i][j-1]+map[i-1][j]+map[i][j+1])/3 : map[i][j]
elseif(j==(dimsize(map,0)-1)&&i!=(dimsize(map,1)-1)&&i!=0)
map[i][j]=map[i][j]==0 ? (map[i][j-1]+map[i-1][j]+map[i+1][j])/3 : map[i][j]
//else
// map[i][j]=map[i][j]==0 ? (map[i-1][j]+map[i][j-1]+map[i+1][j]+map[i][j+1])/4 : map[i][j]
endif
endfor
endfor
// make/d/o/n=(mapsizex/2,mapsizey/2) mapc
// SetScale /P x xMin, 2*delta, mapc
// SetScale /P y yMin, 2*delta, mapc
// mapc=0
// for(i=0;i<dimsize(mapc,0);i++)
// for(j=0;j<dimsize(mapc,1);j++)
// mapc[i][j]+=map[2*i][2*j]
// mapc[i][j]+=map[2*i+1][2*j]
// mapc[i][j]+=map[2*i][2*j+1]
// mapc[i][j]+=map[2*i+1][2*j+1]
// endfor
// endfor
// display
// appendimage mapc
// ModifyImage mapc ctab= {*,*,ColdWarm,0}
End
This does not seem like a question. Maybe you want to post this into the Code Snippets section instead?
July 8, 2024 at 10:49 pm - Permalink
I see. I have never posted in Code Snippets before. Next time I will consider it. Thank you.
July 9, 2024 at 02:13 am - Permalink