
Plot COMSOL multiphysics PINEM simulation

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
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