diff --git a/src/Makefile b/src/Makefile
index 1a92023..abaf845 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -8,7 +8,7 @@ $(info FFTW=$(FFTW))
 FILTER := $(if $(FILTER),$(FILTER),0)
 $(info FILTER=$(FILTER))
 
-# weighting scheme: number-weighted (0), mass-weighted (1), or volume-weighted (2)
+# weighting scheme: number-weighted (0), mass-weighted (1), volume-weighted (2), or emissivity-weighted (3)
 WEIGHT := $(if $(WEIGHT),$(WEIGHT),0)
 $(info WEIGHT=$(WEIGHT))
 
@@ -16,7 +16,7 @@ $(info WEIGHT=$(WEIGHT))
 KERNEL := $(if $(KERNEL),$(KERNEL),0)
 $(info KERNEL=$(KERNEL))
 
-# weighting scheme for the multiscale filter (0: volume-weighted; 1: mass-weighted)
+# weighting scheme for the multiscale filter (0: volume-weighted; 1: mass-weighted, 2: emissivity-weighted)
 WEIGHT_FILTER := $(if $(WEIGHT_FILTER),$(WEIGHT_FILTER),0)
 $(info WEIGHT_FILTER=$(WEIGHT_FILTER))
 
diff --git a/src/filter.f b/src/filter.f
index 0e3719e..7058fb5 100644
--- a/src/filter.f
+++ b/src/filter.f
@@ -1013,6 +1013,10 @@ subroutine multiscale_filter(nx,ny,nz,nl,npatch,pare,
       REAL DENSI_IN0(0:nmax+1,0:nmay+1,0:nmaz+1)
       REAL DENSI_IN1(NAMRX,NAMRY,NAMRZ,NPALEV)
       COMMON /DENSI/ DENSI_IN0,DENSI_IN1
+#elif weight_filter == 2
+      REAL EMISS_IN0(0:nmax+1,0:nmay+1,0:nmaz+1)
+      REAL EMISS_IN1(NAMRX,NAMRY,NAMRZ,NPALEV)
+      COMMON /EMISS/ EMISS_IN0,EMISS_IN1
 #endif 
 
 *     Auxiliary variables
@@ -1045,7 +1049,7 @@ subroutine multiscale_filter(nx,ny,nz,nl,npatch,pare,
 
       lado0 = nx * dx
 
-*     DENS0, DENS1 PROPORTIONAL TO CELLS MASSES or VOLUME!!!
+*     DENS0, DENS1 PROPORTIONAL TO CELLS MASSES, VOLUME, OR EMISSIVITY!!!
       dens0 = 1.0
       do ir=1,nl
        low1=sum(npatch(0:ir-1))+1
@@ -1078,6 +1082,27 @@ subroutine multiscale_filter(nx,ny,nz,nl,npatch,pare,
           dens1(:,:,:,i) = densi_in1(:,:,:,i) * dens1(:,:,:,i) 
         end do
       end do
+#elif weight_filter == 2
+!$omp parallel do shared(nx,ny,nz,emiss_in0,dens0),
+!$omp+            private(i,j,k), 
+!$omp+            default(none)
+      do k=1,nz 
+      do j=1,ny 
+      do i=1,nx 
+        dens0(i,j,k) = emiss_in0(i,j,k) * dens0(i,j,k)
+      end do 
+      end do 
+      end do
+
+      do ir=1,nl 
+        low1=sum(npatch(0:ir-1))+1
+        low2=sum(npatch(0:ir))
+!$omp parallel do shared(low1,low2,emiss_in1,dens1,ir),
+!$omp+            private(i), default(none)
+        do i=low1,low2
+          dens1(:,:,:,i) = emiss_in1(:,:,:,i) * dens1(:,:,:,i) 
+        end do
+      end do
 #endif 
 
       write(*,*) 'in filter'
diff --git a/src/particle_data.f b/src/particle_data.f
index 7f6d1ab..eb484c8 100644
--- a/src/particle_data.f
+++ b/src/particle_data.f
@@ -15,12 +15,18 @@ MODULE PARTICLE_DATA
         REAL ABVC
 #endif
 
-#if weight_scheme == 2
+#if weight_scheme == 2 || weight_filter == 2
         REAL,ALLOCATABLE::VOL(:)
 #else
         ! Dummy variable 
         REAL VOL
-#endif 
+#endif
+
+#if weight_scheme == 3 || weight_filter == 2
+        REAL,ALLOCATABLE::EMISSIVITY(:)
+#else
+        REAL EMISSIVITY
+#endif
 
         INTEGER,ALLOCATABLE::LIHAL(:)
         INTEGER,ALLOCATABLE::LIHAL_IX(:),LIHAL_JY(:),LIHAL_KZ(:)
diff --git a/src/particles.f b/src/particles.f
index 60fd640..ee86da3 100644
--- a/src/particles.f
+++ b/src/particles.f
@@ -1295,7 +1295,6 @@ SUBROUTINE ERROR_PARTICLES(NX,NY,NZ,NL,NPATCH,PATCHNX,PATCHNY,
         UBAS(1:3,1:3,1:3)=U14(IX-1:IX+1,JY-1:JY+1,KZ-1:KZ+1,IPATCH)
         CALL LININT52D_NEW_REAL(AAA,BBB,CCC,RXBAS,RYBAS,RZBAS,UBAS,
      &                          BASZ)
-
        ELSE
         IF (IX.GT.1.AND.IX.LT.NX.AND.
      &      JY.GT.1.AND.JY.LT.NY.AND.
@@ -2315,13 +2314,27 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
       COMMON /DENSI/ L0,L1
 #endif 
 
+#if weight_filter == 2
+      REAL EMISS0(0:NMAX+1,0:NMAY+1,0:NMAZ+1)
+      REAL EMISS1(NAMRX,NAMRY,NAMRZ,NPALEV)
+      COMMON /EMISS/ EMISS0,EMISS1
+#else
+      ! Dummy variables
+      REAL EMISS0, EMISS1
+#endif 
+      
       INTEGER IX,JY,KZ,IR,I,J,K,IPATCH,LOW1,LOW2,CONTA,KNEIGHBOURS
       INTEGER N1,N2,N3,II,JJ,KK,JPATCH,I1,I2,J1,J2,K1,K2,STEP
       INTEGER NPART_TOT,II1,II2,JJ1,JJ2,KK1,KK2,IIP1,JJP1,KKP1
       REAL DXPA,DYPA,DZPA,BASX,BASY,BASZ,BAS,RBAS,BASXX,BASYY,BASZZ
       REAL MEDIOLADO0,PI,MINX,MAXX,MINY,MAXY,MINZ,MAXZ,FRAC_INT,H_KERN
-      REAL*8 BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,BASMASS
+      REAL*8 BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,BASMASS,BAS8EMISS,BAS8_VOL
       REAL,ALLOCATABLE::DIST(:)
+#if weight_filter == 2
+      REAL,ALLOCATABLE::VOL_DIST(:)
+#else
+      REAL VOL_DIST ! dummy variable
+#endif
       INTEGER,ALLOCATABLE::NEIGH(:),LB(:)
       REAL FUIN,U(2,2,2),UW(2,2,2)
       logical do_cell
@@ -2484,10 +2497,12 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
 
 !$OMP PARALLEL DO SHARED(NX,NY,NZ,STEP,I1,I2,J1,J2,K1,K2,RADX,RADY,
 !$OMP+                   RADZ,U2DM,U3DM,U4DM,L0,U2,U3,U4,MASAP,VOL,
+!$OMP+                   EMISSIVITY,EMISS0,
 !$OMP+                   KNEIGHBOURS,DX,VISC0,ABVC,TREE,XTREE,
 !$OMP+                   YTREE,ZTREE,PI,FLAG_MASS),
 !$OMP+            PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8,
-!$OMP+                    BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,DA,SEARCH),
+!$OMP+                    BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,
+!$OMP+                    BAS8EMISS,BAS8_VOL,DA,SEARCH,VOL_DIST),
 !$OMP+            SCHEDULE(DYNAMIC)!, DEFAULT(NONE)
       DO KZ=1,NZ+1,STEP
       DO JY=1,NY+1,STEP
@@ -2500,6 +2515,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
      &     KZ.LT.INT(K2/STEP)*STEP) CYCLE
 
        ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS))
+#if weight_filter == 2
+       ALLOCATE(VOL_DIST(KNEIGHBOURS))
+#endif       
        DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, 
      &       XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)),
      &       ZQUERY=DBLE(RADZ(KZ)), K=KNEIGHBOURS)
@@ -2509,12 +2527,18 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        IF (DIST(KNEIGHBOURS).GT.DX) THEN
         CONTA=KNEIGHBOURS 
        ELSE 
+#if weight_filter == 2
+        DEALLOCATE(VOL_DIST)
+#endif       
         DEALLOCATE(DIST,NEIGH)
         DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, 
      &        XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)),
      &        ZQUERY=DBLE(RADZ(KZ)), RADIUS=DBLE(DX))
         CONTA = DA%SIZE() 
         ALLOCATE(DIST(CONTA), NEIGH(CONTA))
+#if weight_filter == 2
+        ALLOCATE(VOL_DIST(CONTA))
+#endif       
         DIST = DA%V%VALUES 
         NEIGH = DA%I%VALUES
        END IF
@@ -2524,6 +2548,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        L0(IX,JY,KZ)=H_KERN
 
        CALL KERNEL_FUNC(CONTA,CONTA,H_KERN,DIST)
+#if weight_filter == 2
+       DO I=1,CONTA 
+        VOL_DIST(I)=DIST(I)*VOL(NEIGH(I))
+       END DO
+#endif
 #if weight_scheme == 1
        DO I=1,CONTA
         DIST(I)=DIST(I)*MASAP(NEIGH(I))
@@ -2532,6 +2561,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        DO I=1,CONTA 
         DIST(I)=DIST(I)*VOL(NEIGH(I))
        END DO
+#elif weight_scheme == 3
+       DO I=1,CONTA
+          DIST(I)=DIST(I)*EMISSIVITY(NEIGH(I))
+       END DO
 #endif
 
        BAS8=0.D0
@@ -2540,6 +2573,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        BAS8Z=0.D0
        BAS8M=0.D0
        BASMASS=0.D0
+       BAS8EMISS=0.D0
+       BAS8_VOL=0.D0
        DO I=1,CONTA 
         BAS8=BAS8+DIST(I)
         BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I))
@@ -2547,6 +2582,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
         BAS8Z=BAS8Z+DIST(I)*U4DM(NEIGH(I))
 #if use_filter == 1
         BAS8M=BAS8M+DIST(I)*ABVC(NEIGH(I))
+#endif
+#if weight_filter == 2
+        ! emissivity is always volume weighted
+        BAS8_VOL=BAS8_VOL+VOL_DIST(I)
+        BAS8EMISS=BAS8EMISS+VOL_DIST(I)*EMISSIVITY(NEIGH(I))
 #endif
         !BAS8M=MAX(BAS8M,ABVC(NEIGH(I)))
         BASMASS=BASMASS+MASAP(NEIGH(I))
@@ -2558,10 +2598,15 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        VISC0(IX,JY,KZ)=BAS8M/BAS8
 #endif
        !VISC0(IX,JY,KZ)=BAS8M
-
+#if weight_filter == 2
+       EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8_VOL
+#endif
 
        IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3
 
+#if weight_filter == 2
+       DEALLOCATE(VOL_DIST)
+#endif
        DEALLOCATE(DIST, NEIGH)
       END DO 
       END DO 
@@ -2570,7 +2615,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
       ! Fill the blanks by interpolation
 
 !$OMP PARALLEL DO SHARED(NX,NY,NZ,STEP,I1,I2,J1,J2,K1,K2,RADX,RADY,
-!$OMP+                   RADZ,U2,U3,U4,L0,DX,DY,DZ,VISC0),
+!$OMP+                   RADZ,U2,U3,U4,L0,DX,DY,DZ,VISC0,
+!$OMP+                   EMISS0),
 !$OMP+            PRIVATE(IX,JY,KZ,II,JJ,KK,IIP1,JJP1,KKP1,BASX,BASY,
 !$OMP+                    BASZ),
 !$OMP+            DEFAULT(NONE), SCHEDULE(DYNAMIC)
@@ -2639,6 +2685,17 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
      &                 VISC0(IIP1,JJP1,KKP1)* BASX * BASY * BASZ  
 #endif
 
+#if weight_filter == 2
+       EMISS0(IX,JY,KZ)=EMISS0(II,JJ,KK)  *(1-BASX)*(1-BASY)*(1-BASZ) +
+     &                  EMISS0(IIP1,JJ,KK)*  BASX  *(1-BASY)*(1-BASZ) +
+     &                  EMISS0(II,JJP1,KK)*(1-BASX)* BASY *(1-BASZ) +
+     &                  EMISS0(IIP1,JJP1,KK)* BASX * BASY *(1-BASZ) +
+     &                  EMISS0(II,JJ,KKP1)*(1-BASX)*(1-BASY)* BASZ  +
+     &                  EMISS0(IIP1,JJ,KKP1)* BASX *(1-BASY)* BASZ   +
+     &                  EMISS0(II,JJP1,KKP1)*(1-BASX)* BASY * BASZ   +
+     &                  EMISS0(IIP1,JJP1,KKP1)* BASX * BASY * BASZ
+#endif
+
       END DO 
       END DO 
       END DO
@@ -2652,9 +2709,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
 !$OMP+                   JJ2,KK1,KK2,RADX,RADY,RADZ,U2DM,U3DM,
 !$OMP+                   U4DM,L0,U2,U3,U4,KNEIGHBOURS,DX,VISC0,ABVC,
 !$OMP+                   TREE,XTREE,YTREE,ZTREE,FLAG_MASS,
-!$OMP+                   MASAP,VOL,PI),
+!$OMP+                   MASAP,VOL,EMISSIVITY,EMISS0,PI),
 !$OMP+            PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8,
-!$OMP+                    BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,DA,SEARCH),
+!$OMP+                    BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,
+!$OMP+                    BAS8EMISS,BAS8_VOL,DA,SEARCH,VOL_DIST),
 !$OMP+            SCHEDULE(DYNAMIC)!, DEFAULT(NONE)
       DO KZ=K1,K2+1,STEP
       DO JY=J1,J2+1,STEP
@@ -2667,6 +2725,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
      &     KZ.LT.INT(KK2/STEP)*STEP) CYCLE
 
        ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS))
+#if weight_filter == 2
+       ALLOCATE(VOL_DIST(KNEIGHBOURS))
+#endif       
        DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, 
      &       XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)),
      &       ZQUERY=DBLE(RADZ(KZ)), K=KNEIGHBOURS)
@@ -2676,12 +2737,18 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        IF (DIST(KNEIGHBOURS).GT.DX) THEN
         CONTA=KNEIGHBOURS 
        ELSE 
+#if weight_filter == 2
+        DEALLOCATE(VOL_DIST)
+#endif       
         DEALLOCATE(DIST,NEIGH)
         DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, 
      &        XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)),
      &        ZQUERY=DBLE(RADZ(KZ)), RADIUS=DBLE(DX))
         CONTA = DA%SIZE() 
         ALLOCATE(DIST(CONTA), NEIGH(CONTA))
+#if weight_filter == 2
+        ALLOCATE(VOL_DIST(CONTA))
+#endif       
         DIST = DA%V%VALUES 
         NEIGH = DA%I%VALUES
        END IF
@@ -2691,6 +2758,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        L0(IX,JY,KZ)=H_KERN
 
        CALL KERNEL_FUNC(CONTA,CONTA,H_KERN,DIST)
+#if weight_filter == 2
+       DO I=1,CONTA 
+        VOL_DIST(I)=DIST(I)*VOL(NEIGH(I))
+       END DO
+#endif
 #if weight_scheme == 1
       DO I=1,CONTA
        DIST(I)=DIST(I)*MASAP(NEIGH(I))
@@ -2699,6 +2771,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
       DO I=1,CONTA 
        DIST(I)=DIST(I)*VOL(NEIGH(I))
       END DO
+#elif weight_scheme == 3
+      DO I=1,CONTA 
+       DIST(I)=DIST(I)*EMISSIVITY(NEIGH(I))
+      END DO
 #endif
 
        BAS8=0.D0
@@ -2707,6 +2783,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        BAS8Z=0.D0
        BAS8M=0.D0
        BASMASS=0.D0
+       BAS8EMISS=0.D0
+       BAS8_VOL=0.D0
        DO I=1,CONTA 
         BAS8=BAS8+DIST(I)
         BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I))
@@ -2717,6 +2795,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
 #endif
         !BAS8M=MAX(BAS8M,ABVC(NEIGH(I)))
         BASMASS=BASMASS+MASAP(NEIGH(I))
+#if weight_filter == 2
+        BAS8_VOL=BAS8_VOL+VOL_DIST(I)
+        BAS8EMISS=BAS8EMISS+VOL_DIST(I)*EMISSIVITY(NEIGH(I))
+#endif
        END DO
        
        U2(IX,JY,KZ)=BAS8X/BAS8
@@ -2726,9 +2808,15 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        VISC0(IX,JY,KZ)=BAS8M/BAS8
 #endif
        !VISC0(IX,JY,KZ)=BAS8M
-
+#if weight_filter == 2
+       EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8_VOL
+#endif
+       
        IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3
 
+#if weight_filter == 2
+       DEALLOCATE(VOL_DIST)
+#endif
        DEALLOCATE(DIST, NEIGH)
       END DO 
       END DO 
@@ -2738,7 +2826,7 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
 
 !$OMP PARALLEL DO SHARED(NX,NY,NZ,STEP,II1,II2,JJ1,JJ2,KK1,KK2,RADX,
 !$OMP+                   RADY,RADZ,U2,U3,U4,L0,DX,DY,DZ,I1,I2,J1,J2,K1,
-!$OMP+                   K2,VISC0),
+!$OMP+                   K2,VISC0,EMISS0),
 !$OMP+            PRIVATE(IX,JY,KZ,II,JJ,KK,IIP1,JJP1,KKP1,BASX,BASY,
 !$OMP+                    BASZ),
 !$OMP+            DEFAULT(NONE), SCHEDULE(DYNAMIC)
@@ -2807,6 +2895,17 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
      &                 VISC0(IIP1,JJP1,KKP1)* BASX * BASY * BASZ  
 #endif
 
+#if weight_filter == 2
+       EMISS0(IX,JY,KZ)=EMISS0(II,JJ,KK)  *(1-BASX)*(1-BASY)*(1-BASZ) +
+     &                  EMISS0(IIP1,JJ,KK)*  BASX  *(1-BASY)*(1-BASZ) +
+     &                  EMISS0(II,JJP1,KK)*(1-BASX)* BASY *(1-BASZ) +
+     &                  EMISS0(IIP1,JJP1,KK)* BASX * BASY *(1-BASZ) +
+     &                  EMISS0(II,JJ,KKP1)*(1-BASX)*(1-BASY)* BASZ  +
+     &                  EMISS0(IIP1,JJ,KKP1)* BASX *(1-BASY)* BASZ   +
+     &                  EMISS0(II,JJP1,KKP1)*(1-BASX)* BASY * BASZ   +
+     &                  EMISS0(IIP1,JJP1,KKP1)* BASX * BASY * BASZ       
+#endif
+
       END DO 
       END DO
       END DO
@@ -2819,9 +2918,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
 !$OMP PARALLEL DO SHARED(NX,NY,NZ,II1,II2,JJ1,JJ2,KK1,KK2,RADX,RADY,
 !$OMP+                   RADZ,U2DM,U3DM,U4DM,L0,U2,U3,U4,
 !$OMP+                   KNEIGHBOURS,DX,CR0AMR,VISC0,ABVC,
-!$OMP+                   TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL),
+!$OMP+                   TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL,
+!$OMP+                   EMISSIVITY,EMISS0),
 !$OMP+            PRIVATE(IX,JY,KZ,DIST,NEIGH,CONTA,H_KERN,BAS8,
-!$OMP+                    BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,DA,SEARCH),
+!$OMP+                    BAS8X,BAS8Y,BAS8Z,BAS8M,I,BASMASS,
+!$OMP+                    BAS8EMISS,BAS8_VOL,DA,SEARCH,VOL_DIST),
 !$OMP+            SCHEDULE(DYNAMIC)!, DEFAULT(NONE)
       DO KZ=KK1,KK2
       DO JY=JJ1,JJ2
@@ -2829,6 +2930,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
       DO IX=II1,II2
        IF (CR0AMR(IX,JY,KZ).EQ.1) THEN 
         ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS))
+#if weight_filter == 2
+        ALLOCATE(VOL_DIST(KNEIGHBOURS))
+#endif
         DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, 
      &        XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)),
      &        ZQUERY=DBLE(RADZ(KZ)), K=KNEIGHBOURS)
@@ -2838,12 +2942,18 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
         IF (DIST(KNEIGHBOURS).GT.DX) THEN
          CONTA=KNEIGHBOURS 
         ELSE 
+#if weight_filter == 2
+         DEALLOCATE(VOL_DIST)
+#endif       
          DEALLOCATE(DIST,NEIGH)
          DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, 
      &         XQUERY=DBLE(RADX(IX)), YQUERY=DBLE(RADY(JY)),
      &         ZQUERY=DBLE(RADZ(KZ)), RADIUS=DBLE(DX))
          CONTA = DA%SIZE() 
          ALLOCATE(DISt(CONTA), NEIGH(CONTA))
+#if weight_filter == 2
+         ALLOCATE(VOL_DIST(CONTA))
+#endif       
          DIST = DA%V%VALUES 
          NEIGH = DA%I%VALUES
         END IF
@@ -2853,6 +2963,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
         L0(IX,JY,KZ)=H_KERN
 
         CALL KERNEL_FUNC(CONTA,CONTA,H_KERN,DIST)
+#if weight_filter == 2
+       DO I=1,CONTA 
+        VOL_DIST(I)=DIST(I)*VOL(NEIGH(I))
+       END DO
+#endif
 #if weight_scheme == 1
         DO I=1,CONTA
          DIST(I)=DIST(I)*MASAP(NEIGH(I))
@@ -2861,6 +2976,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
         DO I=1,CONTA 
          DIST(I)=DIST(I)*VOL(NEIGH(I))
         END DO
+#elif weight_scheme == 3
+        DO I=1,CONTA 
+         DIST(I)=DIST(I)*EMISSIVITY(NEIGH(I))
+        END DO
 #endif
 
         BAS8=0.D0
@@ -2869,6 +2988,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
         BAS8Z=0.D0
         BAS8M=0.D0
         BASMASS=0.D0
+        BAS8EMISS=0.D0
+        BAS8_VOL=0.D0
         DO I=1,CONTA 
          BAS8=BAS8+DIST(I)
          BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I))
@@ -2879,6 +3000,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
 #endif
          !BAS8M=MAX(BAS8M,ABVC(NEIGH(I)))
          BASMASS=BASMASS+MASAP(NEIGH(I))
+#if weight_filter == 2
+         BAS8_VOL=BAS8_VOL+VOL_DIST(I)
+         BAS8EMISS=BAS8EMISS+VOL_DIST(I)*EMISSIVITY(NEIGH(I))
+#endif
         END DO
       
         U2(IX,JY,KZ)=BAS8X/BAS8
@@ -2888,9 +3013,15 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
         VISC0(IX,JY,KZ)=BAS8M/BAS8
 #endif
         !VISC0(IX,JY,KZ)=BAS8M
-
+#if weight_filter == 2
+        EMISS0(IX,JY,KZ)=BAS8EMISS/BAS8_VOL
+#endif
+        
         IF (FLAG_MASS.EQ.1) L0(IX,JY,KZ)=BASMASS/(4*PI/3)/H_KERN**3
 
+#if weight_filter == 2
+        DEALLOCATE(VOL_DIST)
+#endif
         DEALLOCATE(DIST, NEIGH)
        END IF
       END DO 
@@ -2910,10 +3041,12 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
 !$OMP PARALLEL DO SHARED(LOW1,LOW2,PATCHNX,PATCHNY,PATCHNZ,CR0AMR1,
 !$OMP+                   RX,RY,RZ,U2DM,U3DM,U4DM,L1,U12,U13,U14,
 !$OMP+                   KNEIGHBOURS,DXPA,DX,VISC1,ABVC,SOLAP,
-!$OMP+                   TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL),
+!$OMP+                   TREE,XTREE,YTREE,ZTREE,FLAG_MASS,PI,MASAP,VOL,
+!$OMP+                   EMISSIVITY,EMISS1),
 !$OMP+            PRIVATE(IPATCH,N1,N2,N3,IX,JY,KZ,DIST,NEIGH,DA,
 !$OMP+                    CONTA,H_KERN,BAS8,BAS8X,BAS8Y,BAS8Z,BAS8M,I,
-!$OMP+                    BASMASS,SEARCH,DO_CELL),
+!$OMP+                    BASMASS,BAS8EMISS,BAS8_VOL,SEARCH,DO_CELL,
+!$OMP+                    VOL_DIST),
 !$OMP+            SCHEDULE(DYNAMIC,1)!, DEFAULT(NONE)
        DO IPATCH=LOW1,LOW2 
             !write(*,*) ir,ipatch
@@ -2941,6 +3074,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
           !Q(3)=RZ(KZ,IPATCH)
           
           ALLOCATE(DIST(KNEIGHBOURS), NEIGH(KNEIGHBOURS))
+#if weight_filter == 2
+          ALLOCATE(VOL_DIST(KNEIGHBOURS))
+#endif
           DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, 
      &          XQUERY=DBLE(RX(IX,IPATCH)), YQUERY=DBLE(RY(JY,IPATCH)),
      &          ZQUERY=DBLE(RZ(KZ,IPATCH)), K=KNEIGHBOURS)
@@ -2950,12 +3086,18 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
           IF (DIST(KNEIGHBOURS).GT.DXPA) THEN
            CONTA=KNEIGHBOURS 
           ELSE 
+#if weight_filter == 2
+           DEALLOCATE(VOL_DIST)
+#endif       
            DEALLOCATE(DIST,NEIGH)
            DA = SEARCH%KNEAREST(TREE, XTREE, YTREE, ZTREE, 
      &           XQUERY=DBLE(RX(IX,IPATCH)), YQUERY=DBLE(RY(JY,IPATCH)),
      &           ZQUERY=DBLE(RZ(KZ,IPATCH)), RADIUS=DBLE(DXPA))
            CONTA = DA%SIZE() 
            ALLOCATE(DIST(CONTA), NEIGH(CONTA))
+#if weight_filter == 2
+           ALLOCATE(VOL_DIST(CONTA))
+#endif
            DIST = DA%V%VALUES 
            NEIGH = DA%I%VALUES
           END IF
@@ -2965,6 +3107,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
           !   write(*,*) ipatch,ix,jy,kz,conta,h_kern,dxpa
           !  end if
           CALL KERNEL_FUNC(CONTA,CONTA,H_KERN,DIST)
+#if weight_filter == 2
+          DO I=1,CONTA 
+           VOL_DIST(I)=DIST(I)*VOL(NEIGH(I))
+          END DO
+#endif
 #if weight_scheme == 1
           DO I=1,CONTA
            DIST(I)=DIST(I)*MASAP(NEIGH(I))
@@ -2973,6 +3120,10 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
           DO I=1,CONTA 
             DIST(I)=DIST(I)*VOL(NEIGH(I))
           END DO
+#elif weight_scheme == 3
+          DO I=1,CONTA 
+            DIST(I)=DIST(I)*EMISSIVITY(NEIGH(I))
+          END DO
 #endif
   
           BAS8=0.D0
@@ -2981,6 +3132,8 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
           BAS8Z=0.D0
           BAS8M=0.D0
           BASMASS=0.D0
+          BAS8EMISS=0.D0
+          BAS8_VOL=0.D0
           DO I=1,CONTA 
            BAS8=BAS8+DIST(I)
            BAS8X=BAS8X+DIST(I)*U2DM(NEIGH(I))
@@ -2991,8 +3144,15 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
 #endif
            !BAS8M=MAX(BAS8M,ABVC(NEIGH(I)))
            BASMASS=BASMASS+MASAP(NEIGH(I))
+#if weight_filter == 2
+           BAS8_VOL=BAS8_VOL+VOL_DIST(I)
+           BAS8EMISS=BAS8EMISS+VOL_DIST(I)*EMISSIVITY(NEIGH(I))
+#endif
           END DO
 
+#if weight_filter == 2
+          DEALLOCATE(VOL_DIST)
+#endif
           DEALLOCATE(DIST, NEIGH)
         
           U12(IX,JY,KZ,IPATCH)=BAS8X/BAS8
@@ -3009,6 +3169,9 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
           VISC1(IX,JY,KZ,IPATCH)=BAS8M/BAS8
 #endif
           !VISC1(IX,JY,KZ,IPATCH)=BAS8M
+#if weight_filter == 2
+          EMISS1(IX,JY,KZ,IPATCH)=BAS8EMISS/BAS8_VOL
+#endif
 
           if (flag_mass.eq.0) then
             L1(IX,JY,KZ,IPATCH)=H_KERN
@@ -3044,6 +3207,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
      &    PATCHX,PATCHY,PATCHZ,PATCHRX,PATCHRY,PATCHRZ,
      &    VISC1(1:NAMRX,1:NAMRY,1:NAMRZ,:),NL)
 #endif
+#if weight_filter == 2
+        CALL SYNC_AMR_FILTER(IR,NPATCH,PARE,PATCHNX,PATCHNY,PATCHNZ,
+     &    PATCHX,PATCHY,PATCHZ,PATCHRX,PATCHRY,PATCHRZ,
+     &    EMISS1(1:NAMRX,1:NAMRY,1:NAMRZ,:),NL)
+#endif
 
         LOW1=SUM(NPATCH(0:IR-1))+1
         LOW2=SUM(NPATCH(0:IR))
@@ -3082,6 +3250,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
              u(1:2,1:2,1:2) = VISC1(I:I+1,J:J+1,K:K+1,IPATCH)
              call finer_to_coarser(u,uw,fuin)
              VISC1(II,JJ,KK,JPATCH) = FUIN
+#endif
+#if weight_filter == 2
+             u(1:2,1:2,1:2) = EMISS1(I:I+1,J:J+1,K:K+1,IPATCH)
+             call finer_to_coarser(u,uw,fuin)
+             EMISS1(II,JJ,KK,JPATCH) = FUIN
 #endif
             else
              uw(1:2,1:2,1:2) = 1.
@@ -3106,6 +3279,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
              u(1:2,1:2,1:2) = VISC1(I:I+1,J:J+1,K:K+1,IPATCH)
              call finer_to_coarser(u,uw,fuin)
              VISC0(II,JJ,KK) = FUIN
+#endif
+#if weight_filter == 2
+             u(1:2,1:2,1:2) = EMISS1(I:I+1,J:J+1,K:K+1,IPATCH)
+             call finer_to_coarser(u,uw,fuin)
+             EMISS0(II,JJ,KK) = FUIN
 #endif
             end if
           END DO
@@ -3136,6 +3314,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        WRITE(*,*) 'Mach min,max',BASX,BASY
       END IF
 #endif
+#if weight_filter == 2
+       CALL P_MINMAX_IR(EMISS0,EMISS1,1,1,NX,NY,NZ,NL,PATCHNX,PATCHNY,
+     &                  PATCHNZ,NPATCH,0,BASX,BASY)
+      write(*,*) 'emissivity min,max',BASX,BASY
+#endif
 
 
       DO IR=1,NL
@@ -3162,6 +3345,11 @@ SUBROUTINE INTERPOLATE_VELOCITIES(NX,NY,NZ,NL,NPATCH,PARE,
        ELSE 
         WRITE(*,*) 'Mach min,max',BASX,BASY
        END IF
+#endif
+#if weight_filter == 2
+       CALL P_MINMAX_IR(EMISS0,EMISS1,1,1,NX,NY,NZ,NL,PATCHNX,PATCHNY,
+     &                  PATCHNZ,NPATCH,IR,BASX,BASY)
+       write(*,*) 'emissivity min,max',BASX,BASY
 #endif
       END DO
 
diff --git a/src/reader.f b/src/reader.f
index 9ce4556..6967883 100644
--- a/src/reader.f
+++ b/src/reader.f
@@ -153,8 +153,11 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta,
           deallocate(abvc)
 #endif
 
-#if weight_scheme == 2
+#if weight_scheme == 2 || weight_filter == 2
           deallocate(vol)
+#endif
+#if weight_scheme == 3 || weight_filter == 2
+          deallocate(emissivity)
 #endif
         end if
 
@@ -166,8 +169,11 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta,
        allocate(abvc(parti))
 #endif
 
-#if weight_scheme == 2
-        allocate(vol(parti))
+#if weight_scheme == 2 || weight_filter == 2
+       allocate(vol(parti))
+#endif
+#if weight_scheme == 3 || weight_filter == 2
+       allocate(emissivity(parti))
 #endif
 
        npart(:)=0
@@ -188,11 +194,11 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta,
 ************************************************************************
 ************************************************************************
 
-#if weight_scheme == 2
+#if weight_scheme == 2 || weight_filter == 2
 !$omp parallel do shared(vol, masap, parti), private(i), default(none)
        do i=1,parti
          vol(i)=masap(i)/vol(i)
-       end do
+      end do
 #endif 
        
        npart(0)=low2 !retrocompatibility with general reader
@@ -233,6 +239,10 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta,
         end if
        end if
 #endif
+#if weight_scheme == 3 || weight_filter == 2
+       write(*,*) 'emis=',minval(emissivity(low1:low2)),
+     &      maxval(emissivity(low1:low2))
+#endif
 
        if (xmin.lt.ddxl.or.xmax.gt.ddxr.or.
      &     ymin.lt.ddyl.or.ymax.gt.ddyr.or.
@@ -371,6 +381,10 @@ subroutine read_particles(iter,files_per_snap,nx,ny,nz,t,zeta,
         end if
        end if
 #endif
+#if weight_scheme == 3 || weight_filter == 2
+       write(*,*) 'emis=',minval(emissivity(low1:low2)),
+     &      maxval(emissivity(low1:low2))
+#endif
        
 
        write(*,*) 'routine create mesh ------------------------------'
diff --git a/src/readers/arepo_hdf5.f b/src/readers/arepo_hdf5.f
index cdb0fe9..a828b4b 100644
--- a/src/readers/arepo_hdf5.f
+++ b/src/readers/arepo_hdf5.f
@@ -241,7 +241,7 @@ subroutine read_arepo_hdf5(iter, files_per_snap,
       end if
 #endif
 
-#if weight_scheme == 2
+#if weight_scheme == 2 || weight_filter == 2
       allocate(scr4(NumPart_ThisFile(1)))
       write(*,*) 'reading density ...'
       call h5dopen_f(group_id, "Density", attr_id, status)
@@ -251,6 +251,29 @@ subroutine read_arepo_hdf5(iter, files_per_snap,
       call h5dclose_f(attr_id, status)
       deallocate(scr4)
 #endif
+#if weight_scheme == 3
+      allocate(scr4(numpart_thisfile(1)))
+      write(*,*) 'reading density ...'
+      call h5dopen_f(group_id, "density", attr_id, status)
+      call h5dget_type_f(attr_id, memtype_id, status)
+      call h5dread_f(attr_id, memtype_id, scr4, dims1d, status)
+      emissivity(low1:low2)=scr4(1:numpart_thisfile(1))
+      call h5dclose_f(attr_id, status)
+      write(*,*) 'reading internal energy ...'
+      call h5dopen_f(group_id, "InternalEnergy", attr_id, status)
+      call h5dget_type_f(attr_id, memtype_id, status)
+      call h5dread_f(attr_id, memtype_id, scr4, dims1d, status)
+!$omp parallel do shared(numpart_thisfile, scr4, emissivity, low1), 
+!$omp+            private(i), default(none)
+      do i=1,numpart_thisfile(1)
+         ! calculate the weight as rho^2*sqrt(T)
+         ! as this is only used as weight, we don't convert values to actual temperatures!
+         emissivity(low1+i-1) = emissivity(low1+i-1)*
+     &        emissivity(low1+i-1)*sqrt(scr4(i))
+      end do
+      call h5dclose_f(attr_id, status)
+      deallocate(scr4)
+#endif
         
       call h5gclose_f(group_id, status)
       call h5fclose_f(file_id, status)
@@ -259,4 +282,4 @@ subroutine read_arepo_hdf5(iter, files_per_snap,
 
       return 
       end 
-***********************************************************************
\ No newline at end of file
+***********************************************************************
diff --git a/src/readers/fixgrid_hdf5.f b/src/readers/fixgrid_hdf5.f
index d98a6d2..5eacff3 100644
--- a/src/readers/fixgrid_hdf5.f
+++ b/src/readers/fixgrid_hdf5.f
@@ -109,6 +109,10 @@ subroutine read_fixgrid_hdf5_fields(iter, flag_filter, nx, ny,
        real dens0(0:nmax+1,0:nmay+1,0:nmaz+1)
        real dens1(namrx,namry,namrz,npalev)
        common /densi/ dens0,dens1
+#elif weight_filter == 2
+       real emis0(0:nmax+1,0:nmay+1,0:nmaz+1)
+       real emis1(namrx,namry,namrz,npalev)
+       common /emiss/ emis0,emis1
 #endif 
 
        ! I assume input data is simple precision.
@@ -149,6 +153,16 @@ subroutine read_fixgrid_hdf5_fields(iter, flag_filter, nx, ny,
        end do
        end do
        end do
+#elif weight_filter == 2
+!$omp parallel do shared(emis0,nx,ny,nz), 
+!$omp+ private(ix,jy,kz), default(none)
+       do kz=0,nz+1 
+       do jy=0,ny+1
+       do ix=0,nx+1
+        emis0(ix,jy,kz) = 0.0 
+       end do
+       end do
+       end do
 #endif
 #endif
 
@@ -217,6 +231,31 @@ subroutine read_fixgrid_hdf5_fields(iter, flag_filter, nx, ny,
          call CtoF_order_3d_real_parallel(nx,ny,nz,scr3d)
          dens0(1:nx,1:ny,1:nz) = scr3d(1:nx,1:ny,1:nz)
          call h5dclose_f(attr_id, status)
+#elif weight_filter == 2
+*        Read density and internal energy (if necessary: filter + emissivity-weighting)
+         write(*,*) 'Reading density'
+         call h5dopen_f(file_id, "density", attr_id, status)
+         call h5dget_type_f(attr_id, memtype_id, status)
+         call h5dread_f(attr_id, memtype_id, scr3d, dims3d, status)
+         call CtoF_order_3d_real_parallel(nx,ny,nz,scr3d)
+         emis0(1:nx,1:ny,1:nz) = scr3d(1:nx,1:ny,1:nz)
+         call h5dclose_f(attr_id, status)
+         write(*,*) 'Reading internal energy'
+         call h5dopen_f(file_id, "InternalEnergy", attr_id, status)
+         call h5dget_type_f(attr_id, memtype_id, status)
+         call h5dread_f(attr_id, memtype_id, scr3d, dims3d, status)
+         call CtoF_order_3d_real_parallel(nx,ny,nz,scr3d)
+!$omp parallel do shared(emis0,scr3d,nx,ny,nz), 
+!$omp+ private(ix,jy,kz), default(none)
+         do kz=0,nz+1 
+         do jy=0,ny+1
+         do ix=0,nx+1
+            emis0(ix,jy,kz) = emis(ix,iy,iz)*emis(ix,iy,iz)*
+     &           sqrt(scr3d(ix,iy,iz)) 
+         end do
+         end do
+         end do
+         call h5dclose_f(attr_id, status)
 #endif
 
        end if
@@ -244,6 +283,10 @@ subroutine read_fixgrid_hdf5_fields(iter, flag_filter, nx, ny,
         call p_minmax_ir(dens0,dens1,1,0,nx,ny,nz,nl,patchnx,patchny,
      &                   patchnz,npatch,0,basx,basy)
         write(*,*) 'dens min,max',basx,basy
+#elif weight_filter == 2
+        call p_minmax_ir(emis0,emis1,1,0,nx,ny,nz,nl,patchnx,patchny,
+     &                   patchnz,npatch,0,basx,basy)
+        write(*,*) 'emis min,max',basx,basy
 #endif
         write(*,*)
        end if
@@ -337,4 +380,4 @@ subroutine write_trivial_gridsfile(iter,nx,ny,nz,nl,npatch,
        close(23)
 
        end
-***********************************************************************
\ No newline at end of file
+***********************************************************************
diff --git a/src/readers/gadget_unformatted.f b/src/readers/gadget_unformatted.f
index 38a129c..13dc845 100644
--- a/src/readers/gadget_unformatted.f
+++ b/src/readers/gadget_unformatted.f
@@ -63,7 +63,7 @@ subroutine read_gadget_unformatted(iter, files_per_snap,
 
       integer npart_gadget(6), nall(6), blocksize
       real*8 massarr(6)
-      integer basint1,basint2,basint3,basint4
+      integer basint1,basint2,basint3,basint4,i
       real*8 bas81,bas82,bas83,bas84,bas85,bas86
       integer low1,low2
       real*4,allocatable::scr4(:)
@@ -141,7 +141,7 @@ subroutine read_gadget_unformatted(iter, files_per_snap,
       end if
 #endif
 
-#if weight_scheme == 2
+#if weight_scheme == 2 || weight_filter == 2
       allocate(scr4(sum(npart_gadget(1:6))))
       write(*,*) 'reading density ...'
       call read_float(fil2,'RHO ',scr4,blocksize)
@@ -149,6 +149,25 @@ subroutine read_gadget_unformatted(iter, files_per_snap,
       vol(low1:low2)=scr4(1:npart_gadget(1))
       deallocate(scr4)
 #endif
+#if weight_scheme == 3
+      allocate(scr4(sum(npart_gadget(1:6))))
+      write(*,*) 'reading density ...'
+      call read_float(fil2,'RHO ',scr4,blocksize)
+      write(*,*) ' found for ',(blocksize-8)/4,' particles'
+      emissivity(low1:low2)=scr4(1:npart_gadget(1))
+      write(*,*) 'reading internal energy ...'
+      call read_float(fil2,'U   ',scr4,blocksize)
+      write(*,*) ' found for ',(blocksize-8)/4,' particles'
+!$omp parallel do shared(npart_gadget, scr4, emissivity, low1), 
+!$omp+            private(i), default(none)
+      do i=1,npart_gadget(1)
+         ! calculate the weight as rho^2*sqrt(T)
+         ! as this is only used as weight, we don't convert values to actual temperatures!
+         emissivity(low1+i-1) = emissivity(low1+i-1)*
+     &        emissivity(low1+i-1)*sqrt(scr4(i))
+      end do
+      deallocate(scr4)
+#endif
 
       end do !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/src/readers/masclet.f b/src/readers/masclet.f
index 6141419..1c09f1e 100644
--- a/src/readers/masclet.f
+++ b/src/readers/masclet.f
@@ -148,6 +148,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz,
        real dens0(0:nmax+1,0:nmay+1,0:nmaz+1)
        real dens1(namrx,namry,namrz,npalev)
        common /densi/ dens0,dens1
+#elif weight_filter == 2
+       real emissivity0(0:nmax+1,0:nmay+1,0:nmaz+1)
+       real emissivity1(namrx,namry,namrz,npalev)
+       common /emiss/ emissivity0,emissivity1
 #endif 
 
        integer is_mascletB 
@@ -187,7 +191,17 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz,
         dens0(ix,jy,kz) = 0.0 
        end do
        end do
+      end do
+#elif weight_filter == 2
+!$omp parallel do shared(emissivity0,nx,ny,nz), 
+!$omp+ private(ix,jy,kz), default(none)
+       do kz=0,nz+1 
+       do jy=0,ny+1
+       do ix=0,nx+1
+        emissivity0(ix,jy,kz) = 0.0 
        end do
+       end do
+      end do
 #endif
 #endif
 
@@ -215,6 +229,12 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz,
        do ip = low1,low2 
         dens1(:,:,:,ip) = 0.0 
        end do
+#elif weight_filter == 2
+!$omp parallel do shared(emissivity1,low1,low2)
+!$omp+ private(ip), default(none)
+       do ip = low1,low2 
+        emissivity1(:,:,:,ip) = 0.0 
+       end do
 #endif
 #endif
 
@@ -231,7 +251,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz,
 #if weight_filter == 1
         read(31) (((scr4(i,j,k),i=1,nx),j=1,ny),k=1,nz)
         dens0(1:nx,1:ny,1:nz) = 1.0 + scr4(1:nx,1:ny,1:nz)
-#elif weight_filter == 0 
+#elif weight_filter == 2
+        read(31) (((scr4(i,j,k),i=1,nx),j=1,ny),k=1,nz)
+        emissivity0(1:nx,1:ny,1:nz) = 1.0 + scr4(1:nx,1:ny,1:nz)
+#else
         read(31)
 #endif 
 
@@ -245,7 +268,21 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz,
         read(31) ! pressure
         read(31) ! pot 
         read(31) ! opot 
+#if weight_filter == 2
+        read(31) (((scr4(i,j,k),i=1,nx),j=1,ny),k=1,nz)
+!$omp parallel do shared(emissivity0,scr4,nx,ny,nz), 
+!$omp+ private(ix,jy,kz), default(none)
+        do kz=0,nz+1 
+        do jy=0,ny+1
+        do ix=0,nx+1
+           emissivity0(ix,jy,kz) = emissivity0(ix,jy,kz)*
+     &          emissivity0(ix,jy,kz)*sqrt(scr4(ix,jy,kz)
+        end do
+        end do
+        end do
+#else
         read(31) ! temp 
+#endif
         read(31) ! metal 
         read(31) ! cr0 
         if (is_mascletB.eq.1) then
@@ -267,7 +304,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz,
 #if weight_filter == 1
           read(31) (((scr4(i,j,k),i=1,n1),j=1,n2),k=1,n3)
           dens1(1:n1,1:n2,1:n3,ip) = 1.0 + scr4(1:n1,1:n2,1:n3)
-#elif weight_filter == 0
+#elif weight_filter == 2
+          read(31) (((scr4(i,j,k),i=1,n1),j=1,n2),k=1,n3)
+          emissivity1(1:n1,1:n2,1:n3,ip) = 1.0 + scr4(1:n1,1:n2,1:n3)
+#else
           read(31)
 #endif
           read(31) (((scr4(i,j,k),i=1,n1),j=1,n2),k=1,n3)
@@ -280,6 +320,21 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz,
           read(31) ! pressure
           read(31) ! pot
           read(31) ! opot
+#if weight_filter == 2
+        read(31) (((scr4(i,j,k),i=1,nx),j=1,ny),k=1,nz)
+!$omp parallel do shared(emissivity1,scr4,n1,n2,n3), 
+!$omp+ private(ix,jy,kz), default(none)
+          do kz=0,n1+1
+          do jy=0,n2+1
+          do ix=0,n3+1
+             emissivity1(ix,jy,kz,ip) = emissivity1(ix,jy,kz,ip)*
+     &            emissivity1(ix,jy,kz,ip)*sqrt(scr4(ix,jy,kz)
+          end do
+          end do
+          end do
+#else
+          read(31) ! temp 
+#endif
           read(31) ! temp
           read(31) ! metal
           read(31) ! cr0
@@ -344,6 +399,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz,
       call p_minmax_ir(dens0,dens1,1,0,nx,ny,nz,nl,patchnx,patchny,
      &                 patchnz,npatch,0,basx,basy)
       write(*,*) 'dens min,max',basx,basy
+#elif weight_filter == 2
+      call p_minmax_ir(emis0,emis1,1,0,nx,ny,nz,nl,patchnx,patchny,
+     &                 patchnz,npatch,0,basx,basy)
+      write(*,*) 'emis min,max',basx,basy
 #endif
       call p_minmax_ir(u2,u12,1,1,nx,ny,nz,nl,patchnx,patchny,patchnz,
      &                 npatch,0,basx,basy)
@@ -372,6 +431,10 @@ subroutine read_masclet_fields(iter, flag_filter, nx, ny, nz,
        call p_minmax_ir(dens0,dens1,1,0,nx,ny,nz,nl,patchnx,patchny,
      &                  patchnz,npatch,ir,basx,basy)
        write(*,*) 'dens min,max',basx,basy
+#elif weight_filter == 2
+       call p_minmax_ir(emis0,emis1,1,0,nx,ny,nz,nl,patchnx,patchny,
+     &                  patchnz,npatch,ir,basx,basy)
+       write(*,*) 'emis min,max',basx,basy
 #endif
        call p_minmax_ir(u2,u12,1,1,nx,ny,nz,nl,patchnx,patchny,patchnz,
      &                  npatch,ir,basx,basy)
@@ -450,4 +513,4 @@ subroutine read_masclet_fix_grid_level(nl, nl_vortex, npatch,
        end if
 
        return 
-      end
\ No newline at end of file
+      end