Plot COMSOL multiphysics PINEM simulation

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