diff --git a/csst_ifs_gehong/__init__.py b/csst_ifs_gehong/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..528787cfc8ad81ed41822a8104b60b4896632906 100644
--- a/csst_ifs_gehong/__init__.py
+++ b/csst_ifs_gehong/__init__.py
@@ -0,0 +1 @@
+__version__ = "3.0.0"
diff --git a/csst_ifs_gehong/cube3d.py b/csst_ifs_gehong/cube3d.py
index c7c5c748578974525e383104db5429f35da003a6..c8c7a45266ffd8e7c7a5cba520aa1745459384c2 100644
--- a/csst_ifs_gehong/cube3d.py
+++ b/csst_ifs_gehong/cube3d.py
@@ -2,6 +2,7 @@ import numpy as np
 from .spec1d import StellarContinuum, HII_Region
 import astropy.wcs
 from astropy.io import fits
+from. import __version__
 
 
 class Cube3D():
@@ -31,13 +32,13 @@ class Cube3D():
                                           age=self.stellar_map.age[i, j], feh=self.stellar_map.feh[i, j],
                                           vel=self.stellar_map.vel[i, j], vdisp=self.stellar_map.vdisp[i, j],
                                           ebv=self.stellar_map.ebv[i, j])
+                    self.flux[i, j, :] = self.flux[i, j, :] + ss.flux
                 if self.gas_map is not None:
                     gg = HII_Region(self.config, hii_tem, halpha=self.gas_map.halpha[i, j],
                                     logz=self.gas_map.zh[i, j], vel=self.gas_map.vel[i, j],
                                     vdisp=self.gas_map.vdisp[i, j], ebv=self.gas_map.ebv[i, j])
-                    self.flux[i, j, :] = ss.flux + gg.flux
-                else:
-                    self.flux[i, j, :] = ss.flux
+                    self.flux[i, j, :] = self.flux[i, j, :] + gg.flux
+
 
     def wcs_info(self):
         # wcs = fits.Header()
@@ -69,7 +70,7 @@ class Cube3D():
         hdr = fits.Header()
         hdr['FILETYPE'] = 'SCICUBE'
         hdr['CODE'] = 'csst-ifs-gehong'
-        hdr['VERSION'] = '1.1.3'
+        hdr['VERSION'] = __version__
 
         hdu0 = fits.PrimaryHDU(header=hdr)
 
diff --git a/csst_ifs_gehong/spec1d.py b/csst_ifs_gehong/spec1d.py
index 9b226c02caa4f9e3abb674563153884b23829416..a1535ba3e4cb6263a9cd808eb4dd1566ac588369 100644
--- a/csst_ifs_gehong/spec1d.py
+++ b/csst_ifs_gehong/spec1d.py
@@ -220,7 +220,7 @@ def Calzetti_Law(wave, Rv=4.05):
     reddening_curve = np.zeros(len(wave))
 
     idx = (wave >= 1200) & (wave < 6300)
-    reddening_curve[idx] = 2.659 * (-2.156 + 1.509 * wave_number[idx] - 0.198 * (wave_number[idx] ** 2)) + 0.011 * (wave_number[idx] ** 3) + Rv
+    reddening_curve[idx] = 2.659 * (-2.156 + 1.509 * wave_number[idx] - 0.198 * (wave_number[idx] ** 2) + 0.011 * (wave_number[idx] ** 3)) + Rv
 
     idx = (wave >= 6300) & (wave <= 22000)
     reddening_curve[idx] = 2.659 * (-1.857 + 1.040 * wave_number[idx]) + Rv
@@ -735,13 +735,13 @@ class AGN():
         if (halpha_broad > 0) & (vdisp_broad > 0):
             BLR = AGN_BLR(config, hbeta_flux=halpha_broad / 2.579,
                           hbeta_fwhm=vdisp_broad / 2.355, ebv=ebv, vel=vel)
-            self.flux = self.flux + BLR.flux
+            Fe = AGN_FeII(config, hbeta_broad=halpha_broad / 2.579, ebv=ebv, vel=vel)
+            self.flux = self.flux + BLR.flux + Fe.flux
 
         if (bhmass > 0) & (edd_ratio > 0):
             m5100 = BHmass_to_M5100(bhmass, edd_ratio=edd_ratio, dist=dist)
             PL = AGN_Powerlaw(config, m5100=m5100, ebv=ebv, vel=vel)
-            Fe = AGN_FeII(config, hbeta_broad=halpha_broad / 2.579, ebv=ebv, vel=vel)
-            self.flux = self.flux + PL.flux + Fe.flux
+            self.flux = self.flux + PL.flux
 
 
 def BHmass_to_M5100(bhmass, edd_ratio=0.05, dist=21.0):