init

file_CK04.par

+8
+0.000055 spec/fm25k2odfnew.fits
+0.00017 spec/fm20k2odfnew.fits
+0.00055 spec/fm15k2odfnew.fits
+0.00173 spec/fm10k2odfnew.fits
+0.00546 spec/fm05k2odfnew.fits
+0.01706 spec/fp00k2odfnew.fits
+0.02677 spec/fp02k2odfnew.fits
+0.05202 spec/fp05k2odfnew.fits

main_singlestar.c

+/********************************************************/
+/* code for finding the spec for each of the star  */
+/* within the catalogue                            */
+/* compilation:                                    */
+/* gcc -g -O0 -Wall -I/usr/include -I/usr/include/cfitsio -L/lib64 -lcfitsio -lm main.c -o main */
+/* usage: */
+/* code/main file.par cat.list spec/Ext_odonnell94_R3.1_CK04W.fits */
+/*********************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <fitsio.h>
+#include <time.h>
+#include "utils.c"
+// #include "csst_spec_interp.h"
+
+
+//#include <stdio.h>
+//#include <stdlib.h>
+//#include <math.h>
+#define MAXSTRLEN 200
+//#define G_const 6.67384E-8
+#define lgG log10(6.67384E-8)
+#define lgMsun log10(1.9891E33)
+#define lgExp log10(exp(1.))
+#define twolgpc2cm 2.*log10(3.086e+18)
+
+struct SPECLIB
+{
+  float Z;
+  char name[MAXSTRLEN];
+  int nspec;
+  int nwv;
+  float *logts;
+  float *UniqueLogts;
+  int *iLogts;
+  int NULogts;
+  float *loggs;
+  float *UniqueLoggs;
+  int *iLoggs;
+  int NULoggs;
+  float *wvs;
+  float **specs;
+};
+
+struct STAR
+{
+  float logte;
+  float logg;
+  //float logL;
+  float Mass;
+  float Av;
+  float mu0;
+  float Z;
+  //float mbolmag;
+};
+
+void replace(char * o_string, char * s_string, char * r_string);
+void SortedUniqueElements(float in_array[], float out_array[], int size, int *NU);
+int interp1D (float x, float *v, int *n, int nmax);
+float locateitin1Dtab (float x, float *v, int n, int *i1, int *i2,
+		       int flagrange, int flaglimits, int flaglog);
+
+float locateitin1Dtab(float, float *, int, int *, int *, int, int, int);
+
+void printerror(int);
+
+int loadSpecLibs(char *, char *);
+void loadExts(char *);
+void interpSingleStar(struct STAR, float *, float *);
+
+int NZ=0,NWV=0;
+struct SPECLIB *speclibs=NULL;
+float *Exts=NULL;
+//float *outspec=NULL;
+
+
+void printerror( int status)
+{
+  /*****************************************************/
+  /* Print out cfitsio error messages and exit program */
+  /*****************************************************/
+  if (status)
+    {
+      fits_report_error(stderr, status); /* print error report */
+      exit( status );    /* terminate the program, returning error status */
+    }
+  return;
+}
+
+
+
+
+int loadSpecLibs(char *file_par, char *specDir)
+{
+  fitsfile *fp_spec;
+  FILE *fp_par;
+  char tmpstr[MAXSTRLEN];
+  int iZ,icol,imod,iwv, nspec_table, itable, icol_start, ispec, Nbyte=4;
+  int status=0,nfound,hdutype,nullval,anynulls;
+  long naxes[2];
+  
+  //load the spec libs
+  fp_par=fopen(file_par,"r");
+  fscanf(fp_par, "%d", &NZ);
+  printf("%d different metallicity tables of spec lib\n",NZ);
+  fgets(tmpstr, MAXSTRLEN, fp_par);
+  speclibs=malloc(sizeof(struct SPECLIB)*NZ);
+  for (iZ=0; iZ<NZ; iZ++)
+    {
+      fscanf(fp_par, "%f %s", &speclibs[iZ].Z, speclibs[iZ].name);
+      printf("loading spec file: %s\n",speclibs[iZ].name);
+      fgets(tmpstr, MAXSTRLEN, fp_par);
+      char tmpName[MAXSTRLEN] = "";
+      strcat(tmpName, specDir);
+      strcat(tmpName, speclibs[iZ].name);
+      printf("%s\n",tmpName);
+      if (fits_open_file(&fp_spec, tmpName, READONLY, &status)) printerror(status);
+      if (fits_movabs_hdu(fp_spec, 2, &hdutype, &status)) printerror(status);
+      if (fits_read_keys_lng(fp_spec, "NAXIS", 1, 2, naxes, &nfound, &status)) printerror(status);
+      printf("first table of each metallicity: naxes[0:1]: %ld, %ld\n", naxes[0], naxes[1]);
+      speclibs[iZ].nspec=naxes[1];
+      speclibs[iZ].logts=malloc(sizeof(float)*naxes[1]);
+      speclibs[iZ].UniqueLogts=malloc(sizeof(float)*naxes[1]);
+      speclibs[iZ].iLogts=malloc(sizeof(int)*naxes[1]);
+      speclibs[iZ].loggs=malloc(sizeof(float)*naxes[1]);
+      speclibs[iZ].UniqueLoggs=malloc(sizeof(float)*naxes[1]);
+      speclibs[iZ].iLoggs=malloc(sizeof(int)*naxes[1]);
+      icol=1;
+      fits_read_col(fp_spec, TFLOAT, icol, 1, 1, naxes[1], &nullval, speclibs[iZ].logts, &anynulls, &status);
+      for (imod=0; imod<speclibs[iZ].nspec; imod++) speclibs[iZ].logts[imod]=log10(speclibs[iZ].logts[imod]);
+      fits_read_col(fp_spec, TFLOAT, icol+1, 1, 1, naxes[1], &nullval, speclibs[iZ].loggs, &anynulls, &status);
+
+      SortedUniqueElements(speclibs[iZ].logts, speclibs[iZ].UniqueLogts, speclibs[iZ].nspec, &speclibs[iZ].NULogts);
+      //printf("%d unique Teff grid points of spec lib\n", speclibs[iZ].NULogts);
+      for (imod=0; imod<speclibs[iZ].nspec; imod++)
+	{
+	  for (speclibs[iZ].iLogts[imod]=0;
+	       speclibs[iZ].iLogts[imod]<speclibs[iZ].NULogts;
+	       speclibs[iZ].iLogts[imod]++)
+	    {
+	      if (speclibs[iZ].logts[imod]==speclibs[iZ].UniqueLogts[speclibs[iZ].iLogts[imod]])
+		break;
+	    }
+	}
+      SortedUniqueElements(speclibs[iZ].loggs, speclibs[iZ].UniqueLoggs, speclibs[iZ].nspec, &speclibs[iZ].NULoggs);
+      //printf("%d unique logg grid points of spec lib\n", speclibs[iZ].NULoggs);
+      for (imod=0; imod<speclibs[iZ].nspec; imod++)
+	{
+	  for (speclibs[iZ].iLoggs[imod]=0;
+	       speclibs[iZ].iLoggs[imod]<speclibs[iZ].NULoggs;
+	       speclibs[iZ].iLoggs[imod]++)
+	    {
+	      if (speclibs[iZ].loggs[imod]==speclibs[iZ].UniqueLoggs[speclibs[iZ].iLoggs[imod]])
+		break;
+	    }
+	  //printf("speclibs[iZ].iLoggs[imod]: %d\n",speclibs[iZ].iLoggs[imod]);
+	}
+
+
+
+
+      speclibs[iZ].specs=(float **) malloc(sizeof(float *)*speclibs[iZ].nspec);
+      ispec=0;
+      nspec_table=ceil((speclibs[iZ].nspec+1.)/999.);
+      printf("nspec_table=%d\n", nspec_table);
+      for (itable=0; itable<nspec_table; itable++)
+	{
+	  printf("itable=%d\n", itable);
+	  if (fits_movabs_hdu(fp_spec, 3+itable, &hdutype, &status)) printerror(status);
+	  //if (fits_read_keys_lng(fp_spec, "NAXIS", 1, 2, naxes, &nfound, &status)) printerror(status);
+	  if (fits_read_keys_lng(fp_spec, "NAXIS", 1, 2, naxes, &nfound, &status)) printerror(status);
+	  printf("naxes[0:1]=%ld, %ld\n", naxes[0]/Nbyte, naxes[1]);
+	  icol_start=1;
+	  if (itable==0)
+	    {
+	      speclibs[iZ].nwv=naxes[1];
+	      speclibs[iZ].wvs=malloc(sizeof(float)*naxes[1]);
+	      //fits_read_col(fp_spec, TFLOAT, icol, 1, 1, speclibs[iZ].nwv, &nullval, speclibs[iZ].wvs, &anynulls, &status);
+	      fits_read_col(fp_spec, TFLOAT, 1, 1, 1, speclibs[iZ].nwv, &nullval, speclibs[iZ].wvs, &anynulls, &status);
+	      printf("wavelengths: speclibs[iZ].nwv, speclibs[iZ].wvs[0], speclibs[iZ].wvs[speclibs[iZ].nwv-1]: %d %f %f\n",
+		     speclibs[iZ].nwv, speclibs[iZ].wvs[0], speclibs[iZ].wvs[speclibs[iZ].nwv-1]);
+	      //printf("nspecs: %d\n",speclibs[iZ].nspec);
+	      //if (iZ>1)	{      fits_close_file(fp_spec, &status);      continue;}
+	      icol_start=2;
+	    }
+
+	  //for (icol=icol_start; icol<speclibs[iZ].nspec+2; icol++)
+	  //for (icol=icol_start; icol<naxes[0]/4+1; icol++)
+	  for (icol=icol_start; icol<naxes[0]/Nbyte+1; icol++)
+	    {
+	      //if (ispec==speclibs[iZ].nspec) break;
+	      speclibs[iZ].specs[ispec] = (float *) malloc(sizeof(float)*speclibs[iZ].nwv);
+	      //printf("add: %d",malloc(sizeof(float)*speclibs[iZ].nwv));
+	      fits_read_col(fp_spec, TFLOAT, icol, 1, 1, speclibs[iZ].nwv,
+			    &nullval, speclibs[iZ].specs[ispec], &anynulls, &status);
+	      //printf("%f %f\n",speclibs[iZ].specs[icol-2][0],speclibs[iZ].specs[icol-2][speclibs[iZ].nwv-1]);
+	      for (iwv=0; iwv < speclibs[iZ].nwv; iwv++)
+		{
+		  //speclibs[iZ].specs[icol-2][iwv] = (speclibs[iZ].specs[icol-2][iwv] < 1E-32) ?
+		  //                                   -32 : log10(speclibs[iZ].specs[icol-2][iwv]);
+		  speclibs[iZ].specs[ispec][iwv] = (speclibs[iZ].specs[ispec][iwv] < 1E-20) ?
+		    -20 : log10(speclibs[iZ].specs[ispec][iwv]);
+		}
+	      ispec++;
+	    }
+	  printf("ispec=%d, nspec=%d\n",ispec, speclibs[iZ].nspec);
+	}
+
+
+
+
+
+      
+      fits_close_file(fp_spec, &status);
+      fp_spec=NULL;
+     }
+  NWV=speclibs[0].nwv; //NWV set here and only here
+  //printf("Done loading spec lib files: %d\n",speclibs[0].nwv);
+  return NWV;
+}
+
+
+
+
+void loadExts(char *file_ext)
+{
+  int status=0,hdutype,nullval,anynulls;
+  fitsfile *fp_ext;
+  //float *Exts=NULL;
+  printf("file_ext name is %s\n", file_ext);
+  Exts=malloc(sizeof(float)*NWV);
+  if (fits_open_file(&fp_ext, file_ext, READONLY, &status)) printerror(status);
+  if (fits_movabs_hdu(fp_ext, 2, &hdutype, &status)) printerror(status);
+  fits_read_col(fp_ext, TFLOAT, 2, 1, 1, NWV, &nullval, Exts, &anynulls, &status);
+  printf("Exts[500], Exts[509]: %f %f\n", Exts[500], Exts[509]);
+  fits_close_file(fp_ext, &status);
+  fp_ext=NULL;
+  return;
+}
+
+
+
+
+void interpSingleStar(struct STAR star, float *outspec, float *outspecWave)
+{
+  int iZ,iZsub,iZup,imod;
+  int it1, it2, ig1, ig2;
+  float ftsub, ftup, fgsub, fgup;
+  int iwv;
+  float specsub[NWV], specup[NWV];
+  float wht, totalwht, Zwht;
+  //printf("%f, %f, %f, %f, %f, %f \n",star.Av,star.logg,star.logte, star.Mass, star.mu0, star.Z);
+  //find the two Zs comprise the star
+  iZsub=-99; iZup=-99;
+  if (star.Z<speclibs[0].Z)
+    {
+      iZsub=0; iZup=0;
+    }
+  else if (star.Z>speclibs[NZ-1].Z)
+    {
+      iZsub=NZ-1;iZup=NZ-1;
+    }
+  else
+    {
+      for (iZ=1;iZ<NZ;iZ++)
+	{
+	  if (star.Z<speclibs[iZ].Z)
+	    {
+	      iZsub=iZ-1; iZup=iZ; break;
+	    }
+	}
+    }
+
+  //interp within each Z speclib
+  for (iwv=0; iwv < NWV; iwv++) specsub[iwv]=0.;
+  ftsub = locateitin1Dtab (star.logte, speclibs[iZsub].UniqueLogts,
+			   speclibs[iZsub].NULogts, &it1, &it2, 0, 0, 0);
+  fgsub = locateitin1Dtab (star.logg,  speclibs[iZsub].UniqueLoggs,
+			   speclibs[iZsub].NULoggs, &ig1, &ig2, 0, 0, 0);
+  //printf("sub: it1=%d, it2=%d\n", it1, it2);
+  //printf("sub: ig1=%d, ig2=%d\n", ig1, ig2);
+  totalwht=0.;
+  for (imod=0; imod < speclibs[iZsub].nspec; imod++)
+    {
+      if ( (speclibs[iZsub].iLogts[imod] == it1 || speclibs[iZsub].iLogts[imod] == it2) &&
+	   (speclibs[iZsub].iLoggs[imod] == ig1 || speclibs[iZsub].iLoggs[imod] == ig2) )
+	{
+	  wht = (speclibs[iZsub].iLogts[imod] == it1 ? ftsub : 1.-ftsub) *
+	        (speclibs[iZsub].iLoggs[imod] == ig1 ? fgsub : 1.-fgsub);
+	  totalwht += wht;
+	  for (iwv=0; iwv < speclibs[iZsub].nwv; iwv++)
+	    specsub[iwv] += wht * speclibs[iZsub].specs[imod][iwv];
+	}
+    }
+  //printf("totalwht:%f\n",totalwht);
+  if (totalwht <= 0.) {printf("zero total weight, quit!"); abort();}
+  for (iwv=0; iwv < NWV; iwv++) specsub[iwv] /= totalwht;
+
+  //printf("up-1: it1=%d, it2=%d\n", it1, it2);
+  //printf("up-1: ig1=%d, ig2=%d\n", ig1, ig2);
+  for (iwv=0; iwv < NWV; iwv++) specup[iwv]=0.;
+  ftup = locateitin1Dtab (star.logte, speclibs[iZup].UniqueLogts,
+			  speclibs[iZup].NULogts, &it1, &it2, 0, 0, 0);
+  //printf("up0: it1=%d, it2=%d\n", it1, it2);
+  fgup = locateitin1Dtab (star.logg,  speclibs[iZup].UniqueLoggs,
+			  speclibs[iZup].NULoggs, &ig1, &ig2, 0, 0, 0);
+  //printf("up1: it1=%d, it2=%d\n", it1, it2);
+  //printf("up1: ig1=%d, ig2=%d\n", ig1, ig2);
+  totalwht=0.;
+  for (imod=0; imod < speclibs[iZup].nspec; imod++)
+    {
+      if ( (speclibs[iZup].iLogts[imod] == it1 || speclibs[iZup].iLogts[imod] == it2) &&
+	   (speclibs[iZup].iLoggs[imod] == ig1 || speclibs[iZup].iLoggs[imod] == ig2) )
+	{
+	  wht = (speclibs[iZup].iLogts[imod] == it1 ? ftup : 1.-ftup) *
+	        (speclibs[iZup].iLoggs[imod] == ig1 ? fgup : 1.-fgup);
+	  totalwht += wht;
+	  for (iwv=0; iwv < speclibs[iZup].nwv; iwv++)
+	    specup[iwv] += wht * speclibs[iZup].specs[imod][iwv];
+	}
+    }
+  if (totalwht <= 0.) abort();
+  for (iwv=0; iwv < speclibs[iZup].nwv; iwv++) specup[iwv] /= totalwht;
+
+  Zwht = iZsub == iZup ? 0. : (star.Z - speclibs[iZsub].Z)/(speclibs[iZup].Z - speclibs[iZsub].Z);
+  for (iwv=0; iwv < speclibs[iZsub].nwv; iwv++) //assuming same wv grids
+    outspec[iwv] = specsub[iwv] * Zwht + specup[iwv] * (1.-Zwht);
+
+  //correct for distance and extinction
+  for (iwv=0; iwv < speclibs[iZsub].nwv; iwv++)
+    {
+      outspec[iwv] = outspec[iwv] + lgG + log10(star.Mass) + lgMsun
+	             - star.logg - 0.4*(star.mu0 + 5) - twolgpc2cm;
+      outspec[iwv] = outspec[iwv] - lgExp * Exts[iwv] * star.Av;
+    }
+
+  for (iwv=0; iwv < speclibs[iZsub].nwv; iwv++)
+    {
+      outspecWave[iwv] = speclibs[iZsub].wvs[iwv];
+    }
+  return;
+}
+
+
+

produceSED.py

+'''
+Author: xin zhangxinbjfu@gmail.com
+Date: 2021-08-31 14:58:40
+LastEditors: xin zhangxinbjfu@gmail.com
+LastEditTime: 2022-11-21 16:09:25
+FilePath: /src/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/produceSED_1.py
+Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE
+'''
+
+from tkinter.font import names
+from pylab import *
+import h5py as h5
+from astropy.table import Table
+from scipy import interpolate
+import astropy.constants as cons
+from scipy.interpolate import InterpolatedUnivariateSpline, UnivariateSpline, interp1d
+
+# import healpy as hp
+# from datatable import dt,f
+import numpy as np
+from astropy.cosmology import FlatLambdaCDM
+import os
+
+
+import math
+
+def lyman_forest(wavelen, flux, z):
+    """
+    Compute the Lyman forest mean absorption of an input spectrum,
+    according to D_A and D_B evolution from Madau (1995).
+    The waveeln and flux are in observed frame
+    """
+    if z<=0:
+        flux0 = flux
+    else:
+        nw = 200
+        istep = np.linspace(0,nw-1,nw)
+        w1a, w2a = 1050.0*(1.0+z), 1170.0*(1.0+z)
+        w1b, w2b = 920.0*(1.0+z), 1015.0*(1.0+z)
+        wstepa = (w2a-w1a)/float(nw)
+        wstepb = (w2b-w1b)/float(nw)
+
+        wtempa = w1a + istep*wstepa
+        ptaua = np.exp(-3.6e-03*(wtempa/1216.0)**3.46)
+
+        wtempb = w1b + istep*wstepb
+        ptaub = np.exp(-1.7e-3*(wtempb/1026.0)**3.46\
+                       -1.2e-3*(wtempb/972.50)**3.46\
+                       -9.3e-4*(wtempb/950.00)**3.46)
+
+        da = (1.0/(120.0*(1.0+z)))*np.trapz(ptaua, wtempa)
+        db = (1.0/(95.0*(1.0+z)))*np.trapz(ptaub, wtempb)
+
+        if da>1.0: da=1.0
+        if db>1.0: db=1.0
+        if da<0.0: da=0.0
+        if db<0.0: db=0.0
+        flux0 = flux.copy()
+        id0 = wavelen<=1026.0*(1.0+z)
+        id1 = np.logical_and(wavelen<1216.0*(1.0+z),wavelen>=1026.0*(1.0+z))
+        flux0[id0] = db*flux[id0]
+        flux0[id1] = da*flux[id1]
+
+    return flux0
+def __red(alan,al0,ga,c1,c2,c3,c4):
+
+        fun1 = lambda x: c3/(((x-(al0**2/x))**2)+ga*ga)
+        fun2 = lambda x,cc: cc*(0.539*((x-5.9)**2)+0.0564*((x-5.9)**3))
+        fun = lambda x,cc: c1+c2*x+fun1(x)+fun2(x,cc)
+
+        ala = alan*1.0e-04 #wavelength in microns
+        p = 1.0/ala
+        if np.size(p)>1:
+            p1, p2 = p[p>=5.9], p[p<5.9]
+            ff = np.append(fun(p1,c4), fun(p2,0.0))
+        elif np.size(p)==1:
+            if p<5.9: c4 = 0.0
+            ff = fun(p, c4)
+        else:
+            return
+        return ff
+
+
+
+def reddening(sw, sf, av=0.0, model=0):
+    """
+    calculate the intrinsic extinction of a given template
+
+    Parameters:
+    sw: array
+        wavelength
+    sf: array
+        flux
+    av: float or array
+    model: int
+        Five models will be used:
+        1: Allen (1976) for the Milky Way
+        2: Seaton (1979) fit by Fitzpatrick (1986) for the Milky Way
+        3: Fitzpatrick (1986) for the Large Magellanic Cloud (LMC)
+        4: Prevot et al (1984) and Bouchet (1985) for the Small Magellanic Cloud (SMC)
+        5: Calzetti et al (2000) for starburst galaxies
+        6: Reddy et al (2015) for star forming galaxies
+
+    Return:
+    reddening-corrected flux or observed flux
+    """
+    if model==0 or av==0.0:
+        flux=sf
+    elif model==1: # Allen (1976) for the Milky Way
+        lambda0 = np.array([1000, 1110, 1250, 1430, 1670, \
+                            2000, 2220, 2500, 2850, 3330, \
+                            3650, 4000, 4400, 5000, 5530, \
+                            6700, 9000, 10000, 20000, 100000], dtype=float)
+        kR = np.array([4.20, 3.70, 3.30, 3.00, 2.70, \
+                       2.80, 2.90, 2.30, 1.97, 1.69, \
+                       1.58, 1.45, 1.32, 1.13, 1.00, \
+                       0.74, 0.46, 0.38, 0.11, 0.00],dtype=float)
+        ext0 = InterpolatedUnivariateSpline(lambda0, kR, k=1)
+        A_lambda = av*ext0(sw)
+        A_lambda[A_lambda<0.0] = 0.0
+        flux = sf*10**(-0.4*A_lambda)
+    elif model==2: # Seaton (1979) fit by Fitzpatrick (1986) for the Milky Way
+        Rv=3.1
+        al0, ga, c1, c2, c3, c4 = 4.595, 1.051, -0.38, 0.74, 3.96, 0.26
+        ff11 = __red(1100.0,al0,ga,c1,c2,c3,c4)
+        ff12 = __red(1200.0,al0,ga,c1,c2,c3,c4)
+        slope=(ff12-ff11)/100.0
+        lambda0 = np.array([3650, 4000, 4400, 5000, 5530, \
+                            6700, 9000, 10000, 20000, 100000], dtype=float)
+        kR = np.array([1.58, 1.45, 1.32, 1.13, 1.00, \
+                       0.74, 0.46, 0.38, 0.11, 0.00],dtype=float)
+        fun = interp1d(lambda0, kR, kind='linear')
+
+        sw0 = sw[sw<1200.0]
+        A_lambda0 = (ff11+(sw0-1100.0)*slope)/Rv+1.0
+        sw1 = sw[np.logical_and(sw>=1200.0, sw<=3650.0)]
+        ff = __red(sw1,al0,ga,c1,c2,c3,c4)
+        A_lambda1 = ff/Rv+1.0
+        sw2 = sw[np.logical_and(sw>3650.0, sw<=100000.0)]
+        A_lambda2 = fun(sw2)
+        A_lambda3 = sw[sw>100000.0]*0.0
+        A_lambda = av*np.hstack([A_lambda0,A_lambda1,A_lambda2,A_lambda3])
+        A_lambda[A_lambda<0.0] = 0.0
+        flux = sf*10**(-0.4*A_lambda)
+    elif model==3: # Fitzpatrick (1986) for the Large Magellanic Cloud (LMC)
+        Rv=3.1
+        al0, ga, c1, c2, c3, c4 = 4.608, 0.994, -0.69, 0.89, 2.55, 0.50
+        ff11 = __red(1100.0,al0,ga,c1,c2,c3,c4)
+        ff12 = __red(1200.0,al0,ga,c1,c2,c3,c4)
+        slope=(ff12-ff11)/100.0
+        lambda0 = np.array([3330, 3650, 4000, 4400, 5000, 5530, \
+                            6700, 9000, 10000, 20000, 100000], dtype=float)
+        kR = np.array([1.682, 1.58, 1.45, 1.32, 1.13, 1.00, \
+                       0.74, 0.46, 0.38, 0.11, 0.00],dtype=float)
+        fun = interp1d(lambda0, kR, kind='linear')
+
+        sw0 = sw[sw<1200.0]
+        A_lambda0 = (ff11+(sw0-1100.0)*slope)/Rv+1.0
+        sw1 = sw[np.logical_and(sw>=1200.0, sw<=3330.0)]
+        ff = __red(sw1,al0,ga,c1,c2,c3,c4)
+        A_lambda1 = ff/Rv+1.0
+        sw2 = sw[np.logical_and(sw>3330.0, sw<=100000.0)]
+        A_lambda2 = fun(sw2)
+        A_lambda3 = sw[sw>100000.0]*0.0
+        A_lambda = av*np.hstack([A_lambda0,A_lambda1,A_lambda2,A_lambda3])
+        A_lambda[A_lambda<0.0] = 0.0
+        flux = sf*10**(-0.4*A_lambda)
+    elif model==4: # Prevot et al (1984) and Bouchet (1985) for the Small Magellanic Cloud (SMC)
+        Rv = 2.72
+        lambda0 = np.array([1275, 1330, 1385, 1435, 1490, 1545, \
+                            1595, 1647, 1700, 1755, 1810, 1860, \
+                            1910, 2000, 2115, 2220, 2335, 2445, \
+                            2550, 2665, 2778, 2890, 2995, 3105, \
+                            3704, 4255, 5291, 12500, 16500, 22000], dtype=float)
+        kR = np.array([13.54, 12.52, 11.51, 10.80, 9.84, 9.28, \
+                        9.06, 8.49, 8.01, 7.71, 7.17, 6.90, 6.76, \
+                        6.38, 5.85, 5.30, 4.53, 4.24, 3.91, 3.49, \
+                        3.15, 3.00, 2.65, 2.29, 1.81, 1.00, 0.00, \
+                        -2.02, -2.36, -2.47],dtype=float)
+        kR = kR/Rv+1.0
+        ext0 = InterpolatedUnivariateSpline(lambda0, kR, k=1)
+        A_lambda = av*ext0(sw)
+        A_lambda[A_lambda<0.0] = 0.0
+        flux = sf*10**(-0.4*A_lambda)
+    elif model==5: # Calzetti et al (2000) for starburst galaxies
+        Rv = 4.05
+        sw = sw*1.0e-04           #wavelength in microns
+
+        fun1 = lambda x: 2.659*(-2.156+1.509/x-0.198/x**2+0.011/x**3)+Rv
+        fun2 = lambda x: 2.659*(-1.857+1.040/x)+Rv
+
+        ff11, ff12 = fun1(0.11), fun1(0.12)
+        slope1=(ff12-ff11)/0.01
+        ff99, ff100 = fun2(2.19), fun2(2.2)
+        slope2=(ff100-ff99)/0.01
+
+        sw0 = sw[sw<0.12]
+        sw1 = sw[np.logical_and(sw>=0.12, sw<=0.63)]
+        sw2 = sw[np.logical_and(sw>0.63, sw<=2.2)]
+        sw3 = sw[sw>2.2]
+        k_lambda0 = ff11+(sw0-0.11)*slope1
+        k_lambda1, k_lambda2 = fun1(sw1), fun2(sw2)
+        k_lambda3 = ff99+(sw3-2.19)*slope2
+        A_lambda = av*np.hstack([k_lambda0,k_lambda1,k_lambda2,k_lambda3])/Rv
+        A_lambda[A_lambda<0.0] = 0.0
+        flux = sf*10**(-0.4*A_lambda)
+    elif model==6: # Reddy et al (2015) for satr forming galaxies
+        Rv = 2.505
+        sw = sw*1.0e-04
+
+        fun1 = lambda x: -5.726+4.004/x-0.525/x**2+0.029/x**3+Rv
+        fun2 = lambda x: -2.672-0.010/x+1.532/x**2-0.412/x**3+Rv
+
+        ff11, ff12 = fun1(0.14), fun1(0.15)
+        slope1=(ff12-ff11)/0.01
+        ff99, ff100 = fun2(2.84), fun2(2.85)
+        slope2=(ff100-ff99)/0.01
+
+        sw0 = sw[sw<0.15]
+        sw1 = sw[np.logical_and(sw>=0.15, sw<0.60)]
+        sw2 = sw[np.logical_and(sw>=0.60, sw<2.85)]
+        sw3 = sw[sw>=2.85]
+        k_lambda0 = ff11+(sw0-0.14)*slope1
+        k_lambda1, k_lambda2 = fun1(sw1), fun2(sw2)
+        k_lambda3 = ff99+(sw3-2.84)*slope2
+        A_lambda = av*np.hstack([k_lambda0,k_lambda1,k_lambda2,k_lambda3])/Rv
+        A_lambda[A_lambda<0.0] = 0.0
+        flux  = sf*10**(-0.4*A_lambda)
+
+    else:
+        print("!!! Please select a proper reddening model")
+        sys.exit(0)
+
+    return flux
+
+
+def getObservedSED(sedCat, redshift=0.0, av=0.0, redden=0):
+    z = redshift + 1.0
+    sw, sf = sedCat[:,0], sedCat[:,1]
+    # reddening
+    sf = reddening(sw, sf, av=av, model=redden)
+    sw, sf = sw*z, sf/z
+    # lyman forest correction
+    sf = lyman_forest(sw, sf, redshift)
+    isedObs = (sw.copy(), sf.copy())
+    return isedObs
+
+
+def getSEDData(sedDir='',sedType = 0):
+    sedListF = open(sedDir + 'galaxy.list')
+
+    sedIter = 1
+    l=''
+    while sedIter<=sedType:
+        l = sedListF.readline()
+        sedIter +=1
+    sedfn = l.split()[0]
+
+    sedData = loadtxt(sedDir + sedfn)
+    return sedData
+
+
+def tag_sed(starSpecfile, model_tag, teff=5000, logg=2, feh=0):
+    h5file = h5.File(starSpecfile,'r')
+    model_tag_str = model_tag.decode("utf-8").strip()
+    teff_grid = np.unique(h5file["teff"][model_tag_str])
+    logg_grid = np.unique(h5file["logg"][model_tag_str])
+    feh_grid = np.unique(h5file["feh"][model_tag_str])
+    close_teff = teff_grid[np.argmin(np.abs(teff_grid - teff))]
+    close_logg = logg_grid[np.argmin(np.abs(logg_grid - logg))]
+    if model_tag_str == "koester" or model_tag_str == "MC":
+        close_feh = 99
+    else:
+        close_feh = feh_grid[np.argmin(np.abs(feh_grid - feh))]
+    path = model_tag_str + f"_teff_{close_teff:.1f}_logg_{close_logg:.2f}_feh_{close_feh:.1f}"
+    wave = np.array(h5file["wave"][model_tag_str][()]).ravel()
+    flux = np.array(h5file["sed"][path][()]).ravel()
+    return path, wave, flux
+
+def getABMagAverageVal(ABmag=20.,norm_thr=None, sWave=6840, eWave=8250):
+    """
+    norm_thr: astropy.table, 2 colum, 'WAVELENGTH', 'SENSITIVITY'
+
+    Return:
+        the integerate flux of AB magnitude in the norm_thr(the throughtput of band), photos s-1 m-2 A-1
+    """
+
+    inverseLambda = norm_thr['SENSITIVITY']/norm_thr['WAVELENGTH']
+    norm_thr_i = interpolate.interp1d(norm_thr['WAVELENGTH'], inverseLambda)
+
+    x = np.linspace(sWave,eWave, int(eWave)-int(sWave)+1)
+    y = norm_thr_i(x)
+    AverageLamdaInverse = np.trapz(y,x)/(eWave-sWave)
+    norm = 54798696332.52474 * pow(10.0, -0.4 * ABmag) * AverageLamdaInverse
+    # print('AverageLamdaInverse = ', AverageLamdaInverse)
+    # print('norm = ', norm)
+
+    return norm
+
+def getNormFactorForSpecWithABMAG(ABMag=20., spectrum=None, norm_thr=None, sWave=6840, eWave=8250):
+    """
+    Use AB magnitude system (zero point, fv = 3631 janskys) in the normal band(norm_thr) normalize the spectrum by inpute ABMag
+
+    Parameters
+    ----------
+        spectrum: astropy.table, 2 colum, 'WAVELENGTH', 'FLUX'
+        norm_thr: astropy.table, 2 colum, 'WAVELENGTH', 'SENSITIVITY'
+        sWave: the start of norm_thr
+        eWave: the end of norm_thr
+
+    Return:
+        the normalization factor   flux of AB system(fix a band  and magnitude ) /the flux of inpute spectrum(fix a band) 
+    """
+    spectrumi = interpolate.interp1d(spectrum['WAVELENGTH'], spectrum['FLUX'])
+    norm_thri = interpolate.interp1d(norm_thr['WAVELENGTH'], norm_thr['SENSITIVITY'])
+
+    x = np.linspace(sWave,eWave, int(eWave)-int(sWave)+1)
+
+    y_spec = spectrumi(x)
+    y_thr = norm_thri(x)
+
+    y = y_spec*y_thr
+
+    specAve = np.trapz(y,x)/(eWave-sWave)
+    norm = getABMagAverageVal(ABmag=ABMag, norm_thr=norm_thr, sWave=sWave, eWave=eWave)
+
+    if specAve == 0:
+        return 0
+
+    # print('specAve = ', specAve)
+
+    return norm / specAve
+
+def getABMAG(spec, bandpass_fn, sw=None, ew = None, bpass_type = 'fits'):
+    if bpass_type == 'fits':
+        throughtput = Table.read(bandpass_fn)
+    else:
+        th_dat = np.loadtxt(bandpass_fn)
+        throughtput = Table(th_dat, names=['WAVELENGTH','SENSITIVITY'])
+
+    if sw is None or ew is None:
+        thr_ids = throughtput['SENSITIVITY'] > 0.01
+        sWave = np.floor(throughtput[thr_ids][0][0])
+        eWave = np.ceil(throughtput[thr_ids][-1][0])
+    else:
+        sWave=sw
+        eWave = ew
+        if throughtput['WAVELENGTH'][0]>sWave:
+            throughtput.insert_row(0,[sWave,0.0])
+        if throughtput['WAVELENGTH'][-1]<eWave:
+            throughtput.add_row([eWave,0.0])
+    # print('in getABMAG', sWave, eWave)
+    spectrumi = interpolate.interp1d(spec['WAVELENGTH'], spec['FLUX'])
+    thri = interpolate.interp1d(throughtput['WAVELENGTH'],throughtput['SENSITIVITY'])
+    x = np.linspace(sWave,eWave, (int(eWave)-int(sWave)+1))
+    y_spec = spectrumi(x)*thri(x)
+    interFlux = np.trapz(y_spec, x)
+
+    ABMAG_zero = getABMagAverageVal(
+        ABmag=0, 
+        norm_thr=throughtput, 
+        sWave=sWave, 
+        eWave=eWave)
+    flux_ave = interFlux / (eWave-sWave)
+    ABMAG_spec = -2.5 * np.log10(flux_ave/ABMAG_zero)
+    return ABMAG_spec
+
+def getABMAGANDErr(spec, bandpass_fn, readout = 5, sky = 0.2, dark = 0.02, t = 150, aper = 2, frame = 1, noisepix_num = 22, flux_ratio = 1.0):
+    throughtput = Table.read(bandpass_fn)
+    thr_ids = throughtput['SENSITIVITY'] > 0.01
+    sWave = np.floor(throughtput[thr_ids][0][0])
+    eWave = np.ceil(throughtput[thr_ids][-1][0])
+    # sWave=2000
+    # eWave = 18000
+    # print('in getABMAG', sWave, eWave)
+    spectrumi = interpolate.interp1d(spec['WAVELENGTH'], spec['FLUX'])
+    thri = interpolate.interp1d(throughtput['WAVELENGTH'],throughtput['SENSITIVITY'])
+    x = np.linspace(sWave,eWave, (int(eWave)-int(sWave)+1))
+    y_spec = spectrumi(x)*thri(x)
+    interFlux = np.trapz(y_spec, x)
+
+    ABMAG_zero = getABMagAverageVal(
+        ABmag=0, 
+        norm_thr=throughtput, 
+        sWave=sWave, 
+        eWave=eWave)
+    flux_ave = interFlux / (eWave-sWave)
+    ABMAG_spec = -2.5 * np.log10(flux_ave/ABMAG_zero)
+
+    totalFlux = interFlux * t * frame * math.pi*(aper*aper/4.)
+    noise_var = totalFlux*flux_ratio + (sky + dark) * t * frame * noisepix_num + readout * readout * noisepix_num * frame
+
+    mag_err = 1.087*(noise_var/(totalFlux*flux_ratio))
+
+
+    return ABMAG_spec, mag_err
+
+def getAveEnerge(spec, bandpass_fn):
+    throughtput = Table.read(bandpass_fn)
+    thr_ids = throughtput['SENSITIVITY'] > 0.01
+    sWave = np.floor(throughtput[thr_ids][0][0])
+    eWave = np.ceil(throughtput[thr_ids][-1][0])
+    # sWave=2000
+    # eWave = 18000
+    # print('in getABMAG', sWave, eWave)
+    spectrumi = interpolate.interp1d(spec['WAVELENGTH'], spec['FLUX'])
+    thri = interpolate.interp1d(throughtput['WAVELENGTH'],throughtput['SENSITIVITY'])
+    x = np.linspace(sWave,eWave, (int(eWave)-int(sWave)+1))
+    y_spec = spectrumi(x)
+    interFlux = np.trapz(y_spec, x)
+    return interFlux/(eWave-sWave)
+
+
+def produceSourceSED(sedSoureType = 0,mag_norm = 24.0, norm_filter_thr_fn= 'g.fits',gal_sed_lib_dir = 'Galaxy/', z=0.1, av = 0.1, redden = 0, gal_sedType =1,star_sed_lib_fn='SpecLib.hdf5', lib_tag = 'phoe',teff = 5000, logg = 2, feh = 0): 
+    if sedSoureType == 0:
+        tag = lib_tag.encode('UTF-8')
+        _, wave, flux = tag_sed(star_sed_lib_fn, tag, teff=teff, logg=logg, feh=feh)
+
+    elif sedSoureType==1:
+        sedData = getSEDData(gal_sed_lib_dir, sedType = gal_sedType)
+
+        sed_data = getObservedSED(
+                        sedCat=sedData, 
+                        redshift=z, 
+                        av=av, 
+                        redden=redden)
+        wave = sed_data[0]
+        flux = sed_data[1]
+
+    speci = interpolate.interp1d(wave, flux)
+    lamb = np.arange(2000, 18001 + 0.5, 0.5)
+    y = speci(lamb)
+    # erg/s/cm2/A --> photo/s/m2/A
+    all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+    orig_spec_phot = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+
+    normThr = Table.read(norm_filter_thr_fn)
+    # orig_spec = Table(np.array([wave,flux]).T,  names=(['WAVELENGTH', 'FLUX']))
+    norm_thr_rang_ids = normThr['SENSITIVITY'] > 0.001
+    sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag_norm, spectrum=orig_spec_phot,
+                                                    norm_thr=normThr,
+                                                    sWave=np.floor(normThr[norm_thr_rang_ids][0][0]),
+                                                    eWave=np.ceil(normThr[norm_thr_rang_ids][-1][0]))
+    
+    norm_spec = Table(Table(np.array([wave,flux*sedNormFactor]).T,  names=(['WAVELENGTH', 'FLUX'])))
+    norm_spec_phot = Table(Table(np.array([lamb,all_sed*sedNormFactor]).T,  names=(['WAVELENGTH', 'FLUX'])))
+
+    return norm_spec, norm_spec_phot
+
+def produceNormSED_photon(inputSED = None,mag_norm = 24.0, norm_filter_thr_fn= 'g.fits',ws = 2000, we = 18000):
+    wave = inputSED['WAVELENGTH']
+    flux = inputSED['FLUX']
+    speci = interpolate.interp1d(wave, flux)
+    lamb = np.arange(ws, we + 0.5, 0.5)
+    y = speci(lamb)
+    # erg/s/cm2/A --> photo/s/m2/A
+    all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+    orig_spec_phot = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+
+    normThr = Table.read(norm_filter_thr_fn)
+    # orig_spec = Table(np.array([wave,flux]).T,  names=(['WAVELENGTH', 'FLUX']))
+    norm_thr_rang_ids = normThr['SENSITIVITY'] > 0.001
+    sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag_norm, spectrum=orig_spec_phot,
+                                                    norm_thr=normThr,
+                                                    sWave=np.floor(normThr[norm_thr_rang_ids][0][0]),
+                                                    eWave=np.ceil(normThr[norm_thr_rang_ids][-1][0]))
+    
+    norm_spec_phot = Table(Table(np.array([lamb,all_sed*sedNormFactor]).T,  names=(['WAVELENGTH', 'FLUX'])))
+    return norm_spec_phot
+
+def calculatCSSTMAG_ERR(spec = None, throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/', t = 150,frame = 1, noisepix_num = 22, flux_ratio = 1.0):
+    fil = ['nuv','u','g','r','i','z','y']
+    skybg = {'nuv': 0.00261,'u':0.01823,'g':0.15897,'r':0.20705,'i':0.21433,'z':0.12658,'y':0.037}
+    resMag = {}
+    for fi in fil:
+        mag,err = getABMAGANDErr(spec, throughput_dir+fi+'.Throughput.fits', readout = 5, sky = skybg[fi], dark = 0.02, t = t, aper = 2, frame = frame, noisepix_num = noisepix_num, flux_ratio = flux_ratio)
+        resMag[fi] = [mag,err]
+    return resMag
+
+
+
+def calculatCSSTMAG(spec = None, throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/',ws = None, we = None, filelist = None, band_instr = 'CSST'):
+    fil = ['nuv','u','g','r','i','z','y']
+    if filelist is not None:
+        fil = filelist
+    resMag = {}
+    for fi in fil:
+        if band_instr == 'CSST':
+            mag = getABMAG(spec, throughput_dir+fi+'.Throughput.fits',sw = ws, ew = we, bpass_type = 'fits')
+        else:
+            mag = getABMAG(spec, throughput_dir+fi+'.dat',sw = ws, ew = we, bpass_type = 'dat')
+        resMag[fi] = mag
+    return resMag
+
+def calculatCSSTFilEnergy(spec = None, throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/'):
+    fil = ['nuv','u','g','r','i','z','y']
+    resMag = {}
+    for fi in fil:
+        ene = getAveEnerge(spec, throughput_dir+fi+'.Throughput.fits')
+        resMag[fi] = ene
+    return resMag
+
+def produceGalSED_C6( gal_id_s = '03593100052300144566', gal_z = 1.6927,mag_norm = -99, norm_filter_thr_fn= 'g.fits',galaxy_cat_dir = "/Volumes/EAGET/C6_data/inputData/Catalog_C6_20221212/cat2CSSTSim_bundle/",sedlib_dir = "/Volumes/EAGET/C6_data/inputData/Catalog_C6_20221212/sedlibs/"):
+
+    healPix_id = int(gal_id_s[0:6])
+
+    galcat_file = galaxy_cat_dir + "galaxies_C6_bundle" + gal_id_s[6:12] + '.h5'
+
+    g_id = int(gal_id_s[12:])
+
+    gals_cat = h5.File(galcat_file, 'r')['galaxies']
+
+    coeff = gals_cat[str(healPix_id)]['coeff'][:][g_id]
+
+    pcs = h5.File(os.path.join(sedlib_dir, "pcs.h5"), "r")
+    lamb = h5.File(os.path.join(sedlib_dir, "lamb.h5"), "r")
+    lamb_gal = lamb['lamb'][()]
+    pcs = pcs['pcs'][()]
+
+    cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
+    factor=1
+    if gal_z != 0:
+        factor = 10**(-.4 * cosmo.distmod(gal_z).value)
+    flux = np.matmul(pcs, coeff) * factor
+
+    flux[flux < 0] = 0.
+    sedcat = np.vstack((lamb_gal, flux)).T
+    sed_data = getObservedSED(
+        sedCat=sedcat,
+        redshift=gal_z,
+        av=0.0,
+        redden=0.0
+    )
+    wave, flux = sed_data[0], sed_data[1]
+    speci = interpolate.interp1d(wave, flux)
+    lamb = np.arange(2000, 11001 + 0.5, 0.5)
+    y = speci(lamb)
+    # erg/s/cm2/A --> photo/s/m2/A
+    all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+    orig_spec_phot = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+
+    
+    sedNormFactor=1.0
+    if mag_norm != -99:
+        normThr = Table.read(norm_filter_thr_fn)
+        # orig_spec = Table(np.array([wave,flux]).T,  names=(['WAVELENGTH', 'FLUX']))
+        norm_thr_rang_ids = normThr['SENSITIVITY'] > 0.001
+        sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag_norm, spectrum=orig_spec_phot,
+                                                    norm_thr=normThr,
+                                                    sWave=np.floor(normThr[norm_thr_rang_ids][0][0]),
+                                                    eWave=np.ceil(normThr[norm_thr_rang_ids][-1][0]))
+    norm_spec = Table(Table(np.array([wave,flux*sedNormFactor]).T,  names=(['WAVELENGTH', 'FLUX'])))
+    norm_spec_phot = Table(Table(np.array([lamb,all_sed*sedNormFactor]).T,  names=(['WAVELENGTH', 'FLUX'])))
+
+
+    return norm_spec, norm_spec_phot
+
+# fileDir = os.getcwd()
+
+# #normlization filter star
+# # star_norm_fn = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/SDSS/SLOAN_SDSS.g.fits'
+# star_norm_fn = os.path.join(fileDir, "data/throughputs/SDSS/SLOAN_SDSS.g.fits")
+# #恒星模板库
+# star_sed_lib = "/Volumes/ExtremeSSD/SimData/Catalog_20210126/SpecLib.hdf5"
+
+# #输入参数，星等，得到两个光谱，spec是能量，spec_p是光子
+# spec, spec_photo = produceSourceSED(sedSoureType = 0,mag_norm = 20., norm_filter_thr_fn = star_norm_fn, star_sed_lib_fn=star_sed_lib, lib_tag = 'MM',teff = 3800., logg = 0. , feh = -1.)
+
+
+# #星系星表文件
+# galaxy_cat_dir = "/Volumes/EAGET/C6_data/inputData/Catalog_C6_20221212/cat2CSSTSim_bundle/"
+# #星系光谱模板，PCA
+# sedlib_dir = "/Volumes/EAGET/C6_data/inputData/Catalog_C6_20221212/sedlibs/"
+
+# #normlization filter galaxy
+# # gal_norm_fn = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/LSST/lsst_throuput_g.fits'
+# gal_norm_fn = os.path.join(fileDir, "data/throughputs/LSST/lsst_throuput_g.fits")
+
+# gal_spec, gal_spec_photo = produceGalSED_C6(gal_id_s = '03490800052300010462', gal_z = 0.3764,mag_norm = -99, norm_filter_thr_fn= gal_norm_fn)
+
+
+#根据上面计算的光谱计算csst星等，噪声不需要就不用管了，t, frame, noisepix_num, flux_ratio,都是为了估计噪声
+# mags = calculatCSSTMAG_ERR(spec = spec_photo,throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/', t = 150,frame = 2, noisepix_num = 22, flux_ratio = 1.0)
+
+# # #打印csst星等
+# for k in list(mags.keys()):
+#     print(k, mags[k][0])
+
+# gal_norm_fn = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/LSST/lsst_throuput_g.fits'
+
+# gal_sed_dir = "/Volumes/Extreme SSD/SimData/Templates/Galaxy/"
+
+# # spec1, spec1_p = produceSourceSED(sedSoureType = 1,mag_norm = 22.075, norm_filter_thr_fn= gal_norm_fn,gal_sed_lib_dir = gal_sed_dir, z=0.1, av = 0.1, redden = 0)
+# spec1, spec1_p  = produceSourceSED(sedSoureType = 1,mag_norm = 22.075, norm_filter_thr_fn= gal_norm_fn,gal_sed_lib_dir = gal_sed_dir, z=0.35, av = 0.35, redden = 3.0000,gal_sedType=22)
+# fil = ['nuv','u','g','r','i','z','y']
+# throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/'
+# for fi in fil:
+#     throughput_fn = throughput_dir + fi + '.throughput.fits'
+#     mag = getABMAG(spec_p, throughput_dir+fi+'.Throughput.fits')
+#     print(fi,mag)
+

pycode.py

+'''
+Author: Zhang Xin zhangx@bao.ac.cn
+Date: 2024-01-02 13:34:39
+LastEditors: Zhang Xin zhangx@bao.ac.cn
+LastEditTime: 2024-03-25 13:51:33
+FilePath: /csst-simulation/Users/zhangxin/Work/SlitlessSim/CSST_SIM/Star_spec/csst_spec_interp_clean/code/pycode.py
+Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE
+'''
+##运行下面在mac/linux下执行 
+# cc -fPIC -shared main_singlestar.c -I/usr/include -I/usr/include/cfitsio -lcfitsio -lm -o test.dylib
+
+import os
+import numpy as np
+from astropy.io import fits
+import ctypes
+from astropy.table import Table
+import produceSED
+# from ctypes import *
+
+# struct STAR
+# {
+#   float logte;
+#   float logg;
+#   //float logL;
+#   float Mass;
+#   float Av;
+#   float mu0;
+#   float Z;
+#   //float mbolmag;
+# };
+
+class Star(ctypes.Structure):
+    _fields_ = [
+               ('logte',ctypes.c_float),
+               ('logg',ctypes.c_float),
+               ('Mass',ctypes.c_float),
+               ('Av', ctypes.c_float),
+               ('mu0', ctypes.c_float),
+               ('Z', ctypes.c_float)]
+
+#CHANGE
+#topdir='/run/media/chen/1TS/dupe/gitlab/csst_spec_interp/'
+topdir='/home/zhangxin/CSST_SIM/star_spec/csst_spec_interp_clean/code/'
+code_name='test.so'
+
+d = ctypes.CDLL(os.path.join(topdir,code_name))
+d.loadSpecLibs.argtypes=[ctypes.c_char_p, ctypes.c_char_p]
+d.loadExts.argtypes=[ctypes.c_char_p]
+
+
+#######################################################################################
+#CHANGE
+#nwv = d.loadSpecLibs(str.encode(os.path.join(topdir,'file_CK04.par')),str.encode(topdir))
+#d.loadExts(str.encode(os.path.join(topdir,"spec/Ext_odonnell94_R3.1_CK04W.fits")))
+#TO
+nwv = d.loadSpecLibs(str.encode(os.path.join(topdir,'file_BT-Settl_CSST_wl1000-24000_R1000.par')),str.encode(topdir))
+d.loadExts(str.encode(os.path.join(topdir,"spec/Ext_odonnell94_R3.1_CSST_wl1000-24000_R1000.fits")))
+########################################################################################
+
+print("Done loading Exts")
+spec = (ctypes.c_float*nwv)()
+wave = (ctypes.c_float*nwv)()
+d.interpSingleStar.argtypes=[ctypes.Structure, ctypes.POINTER(ctypes.c_float)]
+
+
+'''
+#example for a single star
+#s=Star(3.620752, 4.701155, 0.599979, 0.067540, 11.200001, 0.008501)
+s=Star(3.9739766, 8.108173, 0.6525841, 0.077022046, 9.05, 0.024376422)
+#s=Star(3.9739766, 4.99, 0.6525841, 0.077022046, 9.05, 0.024376422)
+d.interpSingleStar(s, spec)
+spec_ = spec[:]
+print(spec[500:509])
+'''
+
+from astropy.table import Table
+catalogFn = "/nfsdata/nfsdata/share/CSSTsimInputCat_TH/C9_RA300_DECm60.fits"
+cat = Table.read(catalogFn)
+
+filters = ['nuv','u','g','r','i','z','y']
+filters_other = ['2MASS_H','2MASS_J','2MASS_Ks','GAIA_GAIA3.G','GAIA_GAIA3.Gbp','GAIA_GAIA3.Grp','GALEX_GALEX.FUV','GALEX_GALEX.NUV','LSST_u','LSST_g','LSST_r','LSST_i','LSST_z','LSST_y','PAN-STARRS_PS1.g','PAN-STARRS_PS1.r','PAN-STARRS_PS1.i','PAN-STARRS_PS1.z','PAN-STARRS_PS1.y','SLOAN_SDSS.u','SLOAN_SDSS.g','SLOAN_SDSS.r','SLOAN_SDSS.i','SLOAN_SDSS.z']
+res = {}
+
+nrows = len(cat)
+for fi in filters:
+    res[fi] = np.zeros(nrows)
+for fi in filters_other:
+    res[fi] = np.zeros(nrows)
+
+from mpi4py import MPI
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+rank_size = comm.Get_size()
+print("------------------",rank_size, rank)
+iterNum = 0    
+for star in cat:
+    if iterNum%10000==0:
+        print(iterNum)
+    if iterNum > 1:
+        iterNum = iterNum + 1
+        break
+    if iterNum % rank_size != rank:
+        iterNum = iterNum + 1
+        continue
+    specTable = np.zeros([nwv,2])
+    s=Star(star['mwmsc_logte'], star['mwmsc_logg'], 1., star['mwmsc_av'], star['mwmsc_mu0'], star['mwmsc_z'])
+    #print("star[logTe], star[logg], star[Mass], star[Av], star[mu0], star[Z]: ", star['logTe'], star['logg'], star['Mass'], star['Av'], star['mu0'], star['Z'])
+    d.interpSingleStar(s, spec, wave)
+    specTable[:,0] = wave[:]
+    specTable[:,1] = spec[:]
+    # print(spec[500:509])
+    spec_out = Table(np.array([wave[:], np.power(10,spec[:])]).T, names=('WAVELENGTH', 'FLUX'))
+    spec_norm_phot = produceSED.produceNormSED_photon(inputSED = spec_out,mag_norm = star['mwmsc_gsmag'], norm_filter_thr_fn= 'data/throughputs/SDSS/SLOAN_SDSS.g.fits',ws = 1000, we = 24000)
+    mags = produceSED.calculatCSSTMAG(spec = spec_norm_phot, throughput_dir = 'data/throughputs/CSST/',ws= 2000, we = 11000)
+
+    mags_others = produceSED.calculatCSSTMAG(spec = spec_norm_phot, throughput_dir = 'data/throughputs/filter_transp/',ws= 1000, we = 24000, filelist=filters_other, band_instr='other')
+
+    for fi in filters:
+        res[fi][iterNum] = mags[fi]
+    for fi in filters_other:
+        res[fi][iterNum] = mags_others[fi]
+    iterNum = iterNum + 1
+    print(mags, mags_others)
+    
+
+
+total_res=comm.gather(res, root=0)
+
+if rank == 0:
+    # new_res = res
+    # total_res=comm.gather(res, root=0)
+    # print("gather data size is ", size(total_res))
+    for i in np.arange(1, rank_size, 1):
+        for fi in filters:
+            res[fi] = res[fi] + total_res[i][fi]
+        for fi in filters_other:
+            res[fi] = res[fi] + total_res[i][fi]
+
+    for fi in filters:
+        cat.add_column(np.round(res[fi],5),name='interSpec_'+fi)
+    for fi in filters_other:
+        cat.add_column(np.round(res[fi],5),name='interSpec_'+fi)
+
+    # cat.write("C9_RA300_DECm60_calmag_.fits",overwrite=True)
+    # print('--------------------------')
+    # print(mags['nuv']-star['mwmsc_nuvmag'], mags['u']-star['mwmsc_umag'], mags['g']-star['mwmsc_gmag'], mags['r']-star['mwmsc_rmag'], mags['i']-star['mwmsc_imag'],mags['z']-star['mwmsc_zmag'],mags['y']-star['mwmsc_ymag'])
+#exmple of handling a fits catalouge with many stars
+#wvs=fits.open('spec/fp00k2odfnew.fits')[2].data['wave'] #in Angstrom
+# cat_list=np.loadtxt(topdir+'cat.list', dtype='U')
+# for listfile in cat_list:
+#     hdu = fits.open(topdir+listfile)
+#     cat=hdu[1].data
+#     for star in cat:
+#         specTable = np.zeros([nwv,2])
+#         s=Star(star['logTe'], star['logg'], star['Mass'], star['Av'], star['mu0'], star['Z'])
+#         #print("star[logTe], star[logg], star[Mass], star[Av], star[mu0], star[Z]: ", star['logTe'], star['logg'], star['Mass'], star['Av'], star['mu0'], star['Z'])
+#         d.interpSingleStar(s, spec, wave)
+#         specTable[:,0] = wave[:]
+#         specTable[:,1] = spec[:]
+#         # print(spec[500:509])
+#         spec_out = Table(np.array([wave[:], np.power(10,spec[:])]).T, names=('WAVELENGTH', 'FLUX'))
+#         spec_norm_phot = produceSED.produceNormSED_photon(inputSED = spec_out,mag_norm = 24.0, norm_filter_thr_fn= 'data/throughputs/SDSS/SLOAN_SDSS.g.fits',ws = 2500, we = 10000)
+#         mags = produceSED.calculatCSSTMAG(spec = spec_norm_phot, throughput_dir = 'data/throughputs/CSST/',ws= 2500, we = 10000)
+
+#         # print("MAG is ",mags)
+#         print('--------------------------')
+#         print(mags['nuv']-star['mwmsc_nuvmag'], mag['u']-star['mwmsc_umag'], mag['g']-star['mwmsc_gmag'], mag['r']-star['mwmsc_rmag'], mag['i']-star['mwmsc_imag'],mag['z']-star['mwmsc_zmag'],mag['y']-star['mwmsc_ymag'])
+

readme.md

+cc -fPIC -shared main_singlestar.c -I/home/zhangxin/anaconda3/pkgs/cfitsio-3.470-1/include/ -L/home/zhangxin/anaconda3/pkgs/cfitsio-3.470-1/lib/ -lcfitsio -lm -o test.so

spectres.py

+"""
+SpectRes: A fast spectral resampling function.
+Copyright (C) 2017  A. C. Carnall
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+
+from __future__ import print_function
+import numpy as np 
+import sys
+
+
+# Function for calculating the left hand side (lhs) positions and widths of the spectral bins from their central wavelengths.
+def make_bins(wavelengths, make_rhs="False"):
+    bin_widths = np.zeros(wavelengths.shape[0])
+
+    # This option makes the final entry in the left hand sides array the right hand side of the final bin
+    if make_rhs == "True":
+        bin_lhs = np.zeros(wavelengths.shape[0]+1)
+        #The first lhs position is assumed to be as far from the first central wavelength as the rhs of the first bin.
+        bin_lhs[0] = wavelengths[0] - (wavelengths[1]-wavelengths[0])/2.
+        bin_widths[-1] = (wavelengths[-1] - wavelengths[-2])
+        bin_lhs[-1] = wavelengths[-1] + (wavelengths[-1]-wavelengths[-2])/2.
+        bin_lhs[1:-1] = (wavelengths[1:] + wavelengths[:-1])/2.
+        bin_widths[:-1] = bin_lhs[1:-1]-bin_lhs[:-2]
+
+    # Otherwise just return the lhs positions of each bin
+    else:
+        bin_lhs = np.zeros(wavelengths.shape[0])
+        bin_lhs[0] = wavelengths[0] - (wavelengths[1]-wavelengths[0])/2.
+        bin_widths[-1] = (wavelengths[-1] - wavelengths[-2])
+        bin_lhs[1:] = (wavelengths[1:] + wavelengths[:-1])/2.
+        bin_widths[:-1] = bin_lhs[1:]-bin_lhs[:-1]
+
+    return bin_lhs, bin_widths
+
+
+
+# Function for performing spectral resampling on a spectrum or array of spectra.
+def spectres(spec_wavs, spec_fluxes, resampling, spec_errs=None):
+
+    if (spec_wavs[0]>=resampling[0]):
+        for i, x in enumerate(resampling):
+            if x >= spec_wavs[0]: break
+        if (resampling[0]<1.):
+            print('the output wave min is less than 1A, are you sure????????')
+            print('??????????????????????')
+        if (resampling[0]-(resampling[1]-resampling[0]) > 1.):
+            new_wavs=np.append([resampling[0]-(resampling[1]-resampling[0])],resampling[0:i])
+        else:
+            new_wavs=np.append([1.],resampling[0:i])
+        new_fluxes=0.*new_wavs+1E-100
+        spec_wavs=np.append(new_wavs,spec_wavs)
+        spec_fluxes=np.append(new_fluxes,spec_fluxes)
+
+    if (spec_wavs[-1]<=resampling[-1]):
+        for i, x in enumerate(resampling):
+            if x > spec_wavs[-1]: break
+        new_wavs=np.append(resampling[i:], [resampling[-1]+2.*(resampling[-1]-resampling[-2])])
+        new_fluxes=0.*new_wavs+1E-100
+        spec_wavs=np.append(spec_wavs,new_wavs,)
+        spec_fluxes=np.append(spec_fluxes,new_fluxes)
+
+        
+    # Generate arrays of left hand side positions and widths for the old and new bins
+    filter_lhs, filter_widths = make_bins(resampling, make_rhs="True")
+    spec_lhs, spec_widths = make_bins(spec_wavs)
+
+
+    # Check that the range of wavelengths to be resampled onto falls within the initial sampling region
+    if filter_lhs[0] < spec_lhs[0] or filter_lhs[-1] > spec_lhs[-1]:
+        print("Spec_lhs, filter_lhs, filter_rhs, spec_rhs ", spec_lhs[0], filter_lhs[0], filter_lhs[-1], spec_lhs[-1])
+        #sys.exit("spectres was passed a spectrum which did not cover the full wavelength range of the specified filter curve.")
+        print("spectres was passed a spectrum which did not cover the full wavelength range of the specified filter curve.")
+        
+    
+
+    #Generate output arrays to be populated
+    if spec_fluxes.ndim == 1:
+        resampled = np.zeros((resampling.shape[0]))
+
+    elif spec_fluxes.ndim == 2:
+        resampled = np.zeros((len(resampling), spec_fluxes.shape[1]))
+
+    if spec_errs is not None:
+        if spec_errs.shape != spec_fluxes.shape:
+            sys.exit("If specified, spec_errs must be the same shape as spec_fluxes.")
+        else:
+            resampled_errs = np.copy(resampled)
+
+    start = 0
+    stop = 0
+
+    # Calculate the new spectral flux and uncertainty values, loop over the new bins
+    for j in range(len(filter_lhs)-1):
+
+        # Find the first old bin which is partially covered by the new bin
+        while spec_lhs[start+1] <= filter_lhs[j]:
+            start += 1
+
+        # Find the last old bin which is partially covered by the new bin
+        while spec_lhs[stop+1] < filter_lhs[j+1]:
+            stop += 1
+
+        if spec_fluxes.ndim == 1:
+
+            # If the new bin falls entirely within one old bin the are the same the new flux and new error are the same as for that bin
+            if stop == start:
+
+                resampled[j] = spec_fluxes[start]
+                if spec_errs is not None:
+                    resampled_errs[j] = spec_errs[start]
+
+            # Otherwise multiply the first and last old bin widths by P_ij, all the ones in between have P_ij = 1 
+            else:
+
+                start_factor = (spec_lhs[start+1] - filter_lhs[j])/(spec_lhs[start+1] - spec_lhs[start])
+                end_factor = (filter_lhs[j+1] - spec_lhs[stop])/(spec_lhs[stop+1] - spec_lhs[stop])
+
+                spec_widths[start] *= start_factor
+                spec_widths[stop] *= end_factor
+
+                # Populate the resampled spectrum and uncertainty arrays
+                resampled[j] = np.sum(spec_widths[start:stop+1]*spec_fluxes[start:stop+1])/np.sum(spec_widths[start:stop+1])
+
+                if spec_errs is not None:
+                    resampled_errs[j] = np.sqrt(np.sum((spec_widths[start:stop+1]*spec_errs[start:stop+1])**2.))/np.sum(spec_widths[start:stop+1])
+                
+                # Put back the old bin widths to their initial values for later use
+                spec_widths[start] /= start_factor
+                spec_widths[stop] /= end_factor
+
+
+        # The same as above, except operates on each row of the array, resampling all of the input models
+        elif spec_fluxes.ndim == 2:
+
+            if stop == start:
+
+                resampled[j, :] = spec_fluxes[start, :]
+                if spec_errs is not None:
+                    resampled_errs[j, :] = spec_errs[start, :]
+
+            else:
+
+                start_factor = (spec_lhs[start+1] - filter_lhs[j])/(spec_lhs[start+1] - spec_lhs[start])
+                end_factor = (filter_lhs[j+1] - spec_lhs[stop])/(spec_lhs[stop+1] - spec_lhs[stop])
+
+                spec_widths[start] *= start_factor
+                spec_widths[stop] *= end_factor
+
+                resampled[j, :] = np.sum(np.expand_dims(spec_widths[start:stop+1], axis=1)*spec_fluxes[start:stop+1, :], axis=0)/np.sum(spec_widths[start:stop+1])
+                
+                if spec_errs is not None:
+                    resampled_errs[j, :] = np.sqrt(np.sum((np.expand_dims(spec_widths[start:stop+1], axis=1)*spec_errs[start:stop+1])**2., axis=0))/np.sum(spec_widths[start:stop+1])
+                
+                spec_widths[start] /= start_factor
+                spec_widths[stop] /= end_factor
+
+
+    # If errors were supplied return the resampled spectrum and error arrays
+    if spec_errs is not None:
+        return resampled, resampled_errs
+
+    # Otherwise just return the resampled spectrum array
+    else: 
+        return resampled
+
+

utils.c

+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+/**
+ * The replace function
+ *
+ * Searches all of the occurrences using recursion
+ * and replaces with the given string
+ * @param char * o_string The original string
+ * @param char * s_string The string to search for
+ * @param char * r_string The replace string
+ * @return void The o_string passed is modified
+ */
+void replace(char * o_string, char * s_string, char * r_string) {
+      //a buffer variable to do all replace things
+      char buffer[255];
+      //to store the pointer returned from strstr
+      char * ch;
+
+      //first exit condition
+      if(!(ch = strstr(o_string, s_string)))
+              return;
+
+      //copy all the content to buffer before the first occurrence of the search string
+      strncpy(buffer, o_string, ch-o_string);
+
+      //prepare the buffer for appending by adding a null to the end of it
+      buffer[ch-o_string] = 0;
+
+      //append using sprintf function
+      sprintf(buffer+(ch - o_string), "%s%s", r_string, ch + strlen(s_string));
+
+      //empty o_string for copying
+      o_string[0] = 0;
+      strcpy(o_string, buffer);
+      //pass recursively to replace other occurrences
+      return replace(o_string, s_string, r_string);
+ }
+
+
+
+void SortedUniqueElements(float in_array[], float out_array[], int size, int *NU)
+{
+   int i, j;
+   float temp;
+   *NU=0;
+   //printf("size: %d\n",size);
+   for (i = 0; i < size; i++)
+     {
+       for (j = i+1; j < size; j++) if (in_array[i] == in_array[j]) break;
+       if (j == size) out_array[(*NU)++] = in_array[i];
+     }
+   //printf("NU: %d\n",*NU);
+
+   for ( i = 0; i < (*NU)-1; i++)
+     {
+       for ( j = i+1; j < *NU; j++)
+	 {
+	   if (out_array[i] > out_array[j])
+	     {
+	       temp = out_array[i];
+	       out_array[i] = out_array[j];
+	       out_array[j] = temp;
+	     }
+	 }
+     }
+}
+
+
+
+int  interp1D (float x, float *v, int *n, int nmax)
+{
+  int j, iok; //float tmp;
+  for (j=0, iok=0; j<*n; j++) 
+    if (x==*(v+j)) // if value already exists
+      { iok=1; break; }
+  if (iok==0) // it does not exist, appends value
+    { *(v+j) = x ; (*n)++; }
+  if (*n>nmax)
+    {
+      printf("Error: maximum in dusty table (%d/%d) reached\n", *n, nmax);
+      for (j=0; j<nmax; j++) {printf(" %.3g",*(v+j));} printf("\n");
+      exit(1);
+    }
+  return j; //returns assigned index
+}
+
+
+
+// return interval (i1,i2) and interpolation factor fact on a 1D vector *v
+// aditional parameters: 
+// flagrange= -1 for only negative values, +1 for only positive, 0 for all
+//    (but it requires negative/positive intervals to be contiguous)
+// flaglimits= 0 for no-extrapolation (stay within limits i1,i2), 
+//     -1 for extrapolation in the lower limit
+//     1 for extrapolation in the upper limit
+//     2 for extrapolation in both limits
+float locateitin1Dtab (float x, float *v, int n, int *i1, int *i2,
+		       int flagrange, int flaglimits, int flaglog)
+{
+  int j, lim1, lim2, nlim; 
+  float fact;
+  //linear: flaglog=0; logarithmic: flaglog=1
+  //absolute limits: (lim1,lim2,nlim) instead of (0,n-1,n)
+  lim1=0; lim2=n-1; nlim=n;
+  for (j=0; j<n; j++)
+    {
+      if (flagrange==-1 && *(v+j)>0.0) continue;
+      if (flagrange==+1 && *(v+j)<0.0) continue;
+      lim1=j; break; // first point in interval
+    }
+  for (j=n-1; j>=0; j--)
+    {
+      if (flagrange==-1 && *(v+j)>0.0) continue;
+      if (flagrange==+1 && *(v+j)<0.0) continue;
+      lim2=j; break; // last point in interval
+    }
+  nlim = lim2-lim1;
+  // now search within these limits:
+  if (nlim==1)  // || x <= *(v+lim1)) //PROBLEM HERE?
+    {*i1=lim1; *i2=lim2;}
+  else if  (x == *(v+lim2)) //giada 30/01/2018
+    {*i1=lim2-1; *i2=lim2; }
+  else if (x < *(v+lim1)) //beyond lower limit 
+    {*i1=lim1; *i2=lim1+1;}
+  else if (x > *(v+lim2)) //beyond upper limit 
+    {*i1=lim2-1; *i2=lim2;}
+  else
+    for (j=lim1; j<lim2; j++)
+      {
+	if ( (x >= *(v+j) && x < *(v+j+1)) || //this accepts non-monotonic tables:
+	     (x <= *(v+j) && x > *(v+j+1)) )
+	{ *i1=j; *i2=j+1; break; }
+     }
+
+  if (flaglog == 1)
+     fact = (*i1==*i2) ? 1.0 :
+       (log10(x) - log10(*(v+*i2))) / (log10(*(v+*i1)) - log10(*(v+*i2)));
+  else if (flaglog == 0)
+     fact = (*i1==*i2) ? 1.0 :
+       (x - *(v+*i2)) / (*(v+*i1) - *(v+*i2));
+  else{
+     printf("wrong interpolation flag!\n");
+     exit(1);}
+
+  if (flaglimits==0) // avoids extrapolation
+    {
+      if (fact<0.0) fact=0.0;
+      if (fact>1.0) fact=1.0;
+    }
+  else if (flaglimits==-1) // allows extrapolation to smaller values only 
+    {if (fact>1.0) fact=1.0;}
+  else if (flaglimits==+1) // allows extrapolation to larger values only
+    if (fact<0.0) fact=0.0;
+  // if flaglimits==2, nothing changes -> all extrapolations are allowed
+  //printf ("%f %d %d %f\n", x, *i1, *i2, fact);
+  return fact;
+}