Plot COMSOL multiphysics PINEM simulation

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