Index: /trunk/SophyaLib/SkyT/Makefile
===================================================================
--- /trunk/SophyaLib/SkyT/Makefile	(revision 601)
+++ /trunk/SophyaLib/SkyT/Makefile	(revision 601)
@@ -0,0 +1,40 @@
+# MODULECXXREPNAME := Tests
+# MODULEDECCXXFLAGS := -msg_quiet
+include $(DPCBASEREP)/Include/MakefileUser.h
+# FITSIOLIB2 = -L/exp/planck/FitsIO/cfitsio/OSF1-cxx/
+
+ifdef NOSHLIB
+LIBF =  $(LIB)libSysTools.a $(LIB)libNTools.a $(LIB)libBlitz.a $(LIB)libFitsIOServer.a $(LIB)libPI.a $(LIB)libPIext.a $(LIB)libtools.a 
+LIBS = -L$(LIB) -lSysTools -lNTools -lSamba  -lm
+ifeq ($(CXX),cxx)
+# librairie des instances de templates cxx
+LIBS := $(LIBS) -lertcxx
+endif
+else
+LIBF =  $(SLB)/libtools.so 
+LIBS = -L$(SLB) -ltools -lm 
+ifeq ($(MACHEROS),Linux)
+LIBS := $(LIBS) -ldl
+endif
+endif
+
+# $(OBJ)%.o:%.cc
+
+
+myobjs = $(OBJ)radspec.o $(OBJ)radspecvector.o $(OBJ)nupower.o $(OBJ)blackbody.o $(OBJ)specresp.o $(OBJ)specrespvector.o $(OBJ)squarefilt.o  $(OBJ)trianglefilt.o $(OBJ)gaussfilt.o $(OBJ)convtools.o 
+
+myexe = $(EXE)easyTest $(EXE)skymixer
+
+all : $(myexe)
+
+
+# .PRECIOUS: $(EXE)% $(OBJ)%.o
+$(EXE)easyTest :  $(OBJ)easyTest.o $(myobjs)
+	$(LINK.cc)  $^ $(LIBS) -o $@
+
+$(EXE)skymixer :  $(OBJ)skymixer.o $(myobjs)
+	$(LINK.cc)  $^ $(LIBS) -o $@
+
+
+
+
Index: /trunk/SophyaLib/SkyT/blackbody.cc
===================================================================
--- /trunk/SophyaLib/SkyT/blackbody.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/blackbody.cc	(revision 601)
@@ -0,0 +1,67 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: blackbody.cc,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//      Aim of the class: To give the energy density for a blackbody
+//                        The unity used here is W/m^2/Hz/sr
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+//---------------
+// C++ Headers --
+//---------------
+#include "machdefs.h"
+#include <iostream.h>
+#include <math.h>
+#include "blackbody.h"
+
+//----------------
+// Constructor --
+//----------------
+BlackBody::BlackBody(double temperature)
+  : RadSpectra(0., 10000.)
+{
+  _temperature = temperature;
+}
+
+
+BlackBody::~BlackBody()
+{
+}
+
+
+double
+BlackBody::flux(double nu) const
+{
+  if(nu < -1.e99) nu = -1.e99;
+  if(nu > 1.e99) nu = 1.e99;
+  double temperature = getTemperature();
+  if(nu==0.) return 0.;
+  double hpl = ConvTools::hpl;
+  double cel = ConvTools::cel;
+  double kb  = ConvTools::kb;
+  double puiss1 = nu*pow(10,9);
+  if(puiss1 >  1.e99)   puiss1=1.e99;
+  if(puiss1 < -1.e99)   puiss1=-1.e99;
+  double puiss2 = hpl*nu*pow(10,9)/(kb*temperature);
+  if(puiss2 >  1.e99)   puiss2=1.e99;
+  if(puiss2 < -1.e99)   puiss2=-1.e99;
+  
+  double result=
+   (2*hpl* pow( puiss1 ,3))/( pow(cel,2)*(exp(puiss2)-1));
+  return result;
+}
+
+void
+BlackBody::Print(ostream& os) const
+{
+  os << "BlackBody::Print Temp= " << getTemperature() 
+     << " - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "MeanFreq= " << meanFreq() << "  Emission= " << flux(meanFreq()) << endl;
+  os << "PeakFreq= " << peakFreq() << "  Emission= " << flux(peakFreq()) << endl;
+
+}
Index: /trunk/SophyaLib/SkyT/blackbody.h
===================================================================
--- /trunk/SophyaLib/SkyT/blackbody.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/blackbody.h	(revision 601)
@@ -0,0 +1,43 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: blackbody.h,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+#ifndef BLACKBODY_H_SEEN
+#define BLACKBODY_H_SEEN
+
+#include "radspec.h"
+#include "convtools.h"
+
+class BlackBody : public RadSpectra
+{
+public:   //Constructor
+  //Constructor
+  BlackBody (double temp = ConvTools::tcmb);
+
+  virtual ~BlackBody();
+
+  virtual double flux(double nu) const;
+
+  // returns the value of the temperature
+  inline double getTemperature() const {return _temperature;} 
+
+  // reset the value of the temperature with a new value
+  inline void setTemperature(double newtemp) { _temperature = newtemp; }
+
+  // Redefintion of the Print function
+  virtual void      Print(ostream& os)  const;
+
+protected:
+  double _temperature;
+
+};
+
+#endif
Index: /trunk/SophyaLib/SkyT/convtools.cc
===================================================================
--- /trunk/SophyaLib/SkyT/convtools.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/convtools.cc	(revision 601)
@@ -0,0 +1,24 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: convtools.cc,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+#include "convtools.h"
+
+// Boltzman constant  J/K
+double ConvTools::kb =  1.380662e-23;     
+// Mean value of the cmb temperature (K)
+double ConvTools::tcmb=2.728;
+// constante de Planck en Js
+double ConvTools::hpl=6.626176e-34;   
+// Light speed  PDG 1992 m/s
+double ConvTools::cel=299792458.; 
+
+
+
Index: /trunk/SophyaLib/SkyT/convtools.h
===================================================================
--- /trunk/SophyaLib/SkyT/convtools.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/convtools.h	(revision 601)
@@ -0,0 +1,26 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: convtools.h,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//      Definition of general constant and some convertion functions
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+#ifndef CONV_TOOLS_H_SEEN
+#define CONV_TOOLS_H_SEEN
+
+class ConvTools
+{
+public:
+ static double kb;
+ static double tcmb;
+ static double hpl;
+ static double cel;
+
+};
+#endif
Index: /trunk/SophyaLib/SkyT/detection.cc
===================================================================
--- /trunk/SophyaLib/SkyT/detection.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/detection.cc	(revision 601)
@@ -0,0 +1,334 @@
+//==============================================================================
+/* Nouvelle-Fonction */
+
+// code By Reza Ansari
+
+
+NTuple *  MkAssocNTuple (StarList* stl1, StarList* stl2)
+{
+    // Association des listes d'objets (sources)
+    stl1->Complete();
+    stl2->Complete();
+    StarList::SetPrtLevel(2);
+    stl1->AlGTransfo().SetIdentite(1);
+    int rc = stl1->AlignGeom(stl2, 3, 50, AlG_Transl);
+    cout << " stl1->AlignGeom() -> " << rc << endl;
+    if (rc) { stl1->resTr1_2 = 1;  stl1->AlGTransfo().SetIdentite(1); }
+    int nassoc = stl1->AssocStars(stl2, 3);
+    cout << " stl->AssocStars() NAssoc= " << nassoc << endl;
+    
+    char* nomnt[6] = { "x", "y", "ysz", "yszdet", "dismin", "dism2" } ;
+    NTuple* nt = new NTuple(6, nomnt);
+    float xnt[10];
+    BStar *s1, *s2;
+    int i,j;
+    for(i=0; i<stl1->NbStars(); i++) 
+    {
+	for(j=0; j<10; j++) xnt[j] = -0.;
+	s1 = stl1->Star(i);
+	if (!s1->Nice())  continue;
+	xnt[0] = s1->PosX();
+	xnt[1] = s1->PosY();
+	xnt[2] = s1->Flux();
+	j = stl1->XRef(i);        // probleme ici a cause de aligngeom ou assocstar 
+	if (j >= 0) 
+	{
+	    
+	    s2 = stl2->Star(j);
+	    if (s2->Nice()) 
+	    {
+
+		xnt[3] = s2->Flux();
+		xnt[4] = stl1->DAss(i);
+		xnt[5] = stl1->DAss2(i);
+	    }
+        }
+	nt->Fill(xnt); 
+    }
+    return(nt);
+}
+
+//==============================================================================
+/* Nouvelle-Fonction */ 
+// la fonction repond a la question : le point de coordonnees (x,y) est-
+// il un maximum local ( sur un pave 9*9 centre sur (x,y) )
+// ceci quelques soient x et y
+
+void MaxImg(ImageR4* im,double x,double y,bool& max)
+{
+  
+   	
+    if ( ( x>0 ) && ( x<im->XSize()-1 ) && ( y>0 ) && ( y<im->YSize()-1 ) )
+    // pour le centre de l'image (8 points a verifier)
+    {
+		
+	if (( (*im)(x,y) > (*im)(x-1,y-1) )&& ( (*im)(x,y) > (*im)(x,y-1) ))
+	{
+	    if (( (*im)(x,y) > (*im)(x+1,y-1) ) && ( (*im)(x,y) > (*im)(x-1,y+1) ))
+	    {
+		if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
+		{
+		    if (( (*im)(x,y) > (*im)(x-1,y) ) && ( (*im)(x,y) > (*im)(x+1,y) ))
+		    {
+			max= true;
+		    }
+		}
+	    }
+	}
+    }
+    if ( (x==0) &&  ( y>0 ) && ( y<im->YSize()-1 ) )
+    // bord en x a "gauche" (5 points a verifier)
+    {
+	if (( (*im)(x,y) > (*im)(x+1,y-1) ) && ( (*im)(x,y) > (*im)(x,y-1) ))
+	{
+	    if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
+	    {	
+		if ( (*im)(x,y) > (*im)(x+1,y) )
+		{
+		    max= true;
+		}
+	    }
+	}
+    }
+    if ( (y==0) &&  ( x>0 ) && ( x<im->XSize()-1 ) )
+    // bord en y en haut (5 points )
+    {
+	if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
+	{
+	    if (( (*im)(x,y) > (*im)(x-1,y) ) && ( (*im)(x,y) > (*im)(x+1,y) ))
+	    {
+		if ( (*im)(x,y) > (*im)(x-1,y+1) )
+		{
+		    max= true;
+		}
+	    }
+	}
+    }
+    if ( (x==im->XSize()-1 ) &&  ( y>0 ) && ( y<im->YSize()-1 ) )
+    // bord en x a "droite" (5 points a verifier)
+    {
+	if (( (*im)(x,y) > (*im)(x-1,y-1) ) && ( (*im)(x,y) > (*im)(x,y-1) ))
+	{
+	    if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x-1,y+1) ))
+	    {	
+		if ( (*im)(x,y) > (*im)(x-1,y) )
+	       	{
+		    max= true;
+		}
+	    }
+	}
+    }
+    if ( (y==im->YSize()-1 ) &&  ( x>0 ) && ( x<im->XSize()-1 ) )
+    // bord en y en bas (5 points )
+    {
+	if (( (*im)(x,y) > (*im)(x,y-1) ) && ( (*im)(x,y) > (*im)(x+1,y-1) ))
+	{
+	    if (( (*im)(x,y) > (*im)(x-1,y) ) && ( (*im)(x,y) > (*im)(x+1,y) ))
+	    {
+		if ( (*im)(x,y) > (*im)(x-1,y-1) )
+		{
+		    max= true;
+		}
+	    }
+	}
+    }
+    
+    if ( (x==0) &&  ( y==0 ) )
+    // coin en haut a "gauche" (3 points a verifier)
+    {
+	if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
+	{	
+	    if ( (*im)(x,y) > (*im)(x+1,y) )
+	    {
+		max= true;
+	    }
+	}
+    }
+    
+
+    if ( (x==im->XSize()-1 ) &&  ( y==0 ) )
+    // coin en haut a "droite" (3 points a verifier)
+    {
+	if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x-1,y+1) ))
+	{	
+	    if ( (*im)(x,y) > (*im)(x-1,y) )
+	    {
+		max= true;
+	    }
+	}
+    }
+		
+    if ( (x==im->XSize()-1 ) &&  ( y==im->YSize()-1 ) )
+    // coin en bas a "droite" (3 points a verifier)
+    {
+	if (( (*im)(x,y) > (*im)(x,y-1) ) && ( (*im)(x,y) > (*im)(x-1,y-1) ))
+	{	
+	    if ( (*im)(x,y) > (*im)(x-1,y) )
+	    {
+		max= true;
+	    }
+	}
+    }
+    if ( (x==0) &&  ( y==im->YSize()-1 ) )
+    // coin en bas a "gauche" (3 points a verifier)
+    {
+	if (( (*im)(x,y) > (*im)(x,y-1) ) && ( (*im)(x,y) > (*im)(x+1,y-1) ))
+	{	
+	    if ( (*im)(x,y) > (*im)(x+1,y) )
+	    {
+		max= true;
+	    }
+	}
+    }
+}
+
+
+//==============================================================================
+/* Nouvelle-Fonction */ 
+// calcule la somme des pixels dans un carre 9*9 centre sur (x,y)
+// dans l'image im, ceci quelques soient x et y
+
+void SomPix(ImageR4& im,double x,double y,float& som)
+{
+  
+   	
+    if ( ( x>0 ) && ( x<im.XSize()-1 ) && ( y>0 ) && ( y<im.YSize()-1 ))
+    // pour le centre de l'image (9 points)
+    {
+		
+	som=  im(x,y) +  im(x,y+1) + 
+		im(x,y-1) + im(x+1,y) +
+		im(x-1,y) + im(x-1,y-1) +
+		im(x-1,y+1) + im(x+1,y+1) +
+		im(x+1,y-1);
+    }
+    if ( (x==0) &&  ( y>0 ) && ( y<im.YSize()-1 ) )
+    // bord en x a "gauche" (6 points )
+    {
+	som= im(x,y) +  im(x,y+1) + im(x,y-1) + im(x+1,y) + im(x+1,y+1) +
+		im(x+1,y-1);
+    }
+
+    if ( (y==0) &&  ( x>0 ) && ( x<im.XSize()-1 ) )
+    // bord en y en haut (6 points )
+    {
+	som=  im(x,y) +  im(x,y+1) + im(x+1,y) +
+		im(x-1,y) +	im(x-1,y+1) + im(x+1,y+1) ;
+    }
+
+    if ( (x==im.XSize()-1 ) &&  ( y>0 ) && ( y<im.YSize()-1 ) )
+    // bord en x a "droite" (6 points)
+    {
+	som= im(x,y) +  im(x,y+1) + 
+	       im(x,y-1) +	im(x-1,y) + im(x-1,y-1) +
+	       im(x-1,y+1) ;
+
+    }
+    if ( (y==im.YSize()-1 ) &&  ( x>0 ) && ( x<im.XSize()-1 ) )
+    // bord en y en bas (6 points )
+    {
+	som=  im(x,y) +	im(x,y-1) + im(x+1,y) +
+		im(x-1,y) + im(x-1,y-1) +	im(x+1,y-1);
+    }
+    
+    if ( (x==0) &&  ( y==0 ) )
+    // coin en haut a "gauche" (4 points )
+    {
+	som= im(x,y) +im(x+1,y) + im(x,y+1) + im(x+1,y+1) ;
+    }
+    
+
+    if ( (x==im.XSize()-1 ) &&  ( y==0 ) )
+    // coin en haut a "droite" (4 points )
+    {
+	som= im(x,y) + im(x-1,y) + im(x-1,y+1) + im(x,y+1) ;
+    }
+		
+    if ( (x==im.XSize()-1 ) &&  ( y==im.YSize()-1 ) )
+    // coin en bas a "droite" (4 points )
+    {
+	som= im(x,y) + im(x-1,y) + im(x-1,y-1) + im(x,y-1) ;
+    }
+    if ( (x==0) &&  ( y==im.YSize()-1 ) )
+    // coin en bas a "gauche" (4 points)
+    {
+	som= im(x,y) + im(x+1,y) + im(x+1,y-1) + im(x,y-1) ;
+    }
+    
+
+}
+
+/* ------Nouvelle fonction -----------*/
+// detection de l'effet SZ
+// la carte des sources detectees est mise dans nb_det (d'ou la tramsmission par reference)
+// le nombre de sources detectees est mis dans nb (d'ou la tramsmission par reference)
+// la fonction retourne une StarList
+
+StarList * detec_SZ(ImageR4 & corin, ImageR4& nb_det,double& nb,int taille_f)
+{
+    ImageR4 * cor=NULL; 
+    // On cree une liste d'etoiles pour stocker les sources detectees
+    StarList* stl = new StarList(NULL, 1000, 250); 
+    BStar* star;
+    
+    if (taille_f <= 1 ) cor = &corin;             // si l'on ne veut pas de filtrage il suffit d'appeler
+                                                  // la fonction detec_SZ(..,.,.,1)
+    else 
+    { 
+	ImageR4 filt_c(taille_f,taille_f);        // declaration et ...
+	filt_c=1.;                                // initialisation du filtre 
+	cor = FilterImage(&corin,cor,&filt_c);    // filtrage
+    }
+
+    nb=0;
+    nb_det=0;
+    
+    
+    double seuil=9/taille_f;                        // on garde 3 sigma pour chaque dimension x et y
+    //    cout << "in seuil!!!" << seuil << endl;
+    double taille=cor->XSize()*cor->YSize();
+
+    bool tmp=false;
+    bool& fgdet=tmp;
+    float posx, posy, flux, fond, sompix, pixmax;   // parametre d'un etoile (source SZ)
+    
+
+
+    for (int i=0;i<cor->XSize();i++)    
+    {
+	for (int j=0;j<cor->YSize();j++)
+	{
+            fgdet = false;
+	    
+	    if ( (*cor)(i,j)>seuil )        // image filtree ou non filtree  
+	    { 
+		MaxImg(cor,i,j,fgdet);
+		if ( fgdet )
+		{
+		    
+		    nb=nb+1;          // calcul du nombre de sources
+		    nb_det(i,j)=10;   // valeur arbitraire qui pourra etre remplacee par le y mesure par exemple
+		    //fgdet = true;
+		 
+		}
+	    }
+
+	   
+            if ( !fgdet ) continue;                     // creation de la liste d'etoiles
+            posx = i+0.5;  posy = j+0.5;                      // on utilise les caracteristiques de la carte
+	                                              // non filtree
+            pixmax = corin(i,j);
+	    SomPix(corin,i,j,sompix);
+	    flux = sompix;
+            fond = 0.;
+	    star = stl->Next();
+	    star->SetX(posx, posy, flux, fond, pixmax, sompix, 
+                       posx, posy, kCourRef);
+	    star->ClearFlag(BStar::flagAll, kCourRef);
+	    star->SetFlag(BStar::flagOK, kCourRef);
+	}
+    }
+    return(stl);
+}
+
+
Index: /trunk/SophyaLib/SkyT/easyTest.cc
===================================================================
--- /trunk/SophyaLib/SkyT/easyTest.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/easyTest.cc	(revision 601)
@@ -0,0 +1,135 @@
+#include <iostream.h>
+#include "blackbody.h"
+#include "radspecvector.h"
+#include "specrespvector.h"
+#include "tvector.h"
+// #include "szsource.h"
+#include "squarefilt.h"
+#include "trianglefilt.h"
+int
+main() 
+{ 
+  cout << "BlackBody!" << endl;
+  BlackBody myBB;
+  double freqmin =  10.;
+  double freqmax = 100;
+  double logIF = myBB.logIntegratedFlux(freqmin,freqmax);
+  double IF    = myBB.integratedFlux(freqmin,freqmax);
+  cout << "log   " << logIF << endl;
+  cout << "nolog " << IF    << endl;
+  cout << "min et max " << myBB.minFreq() << " :: " << myBB.maxFreq() << endl;
+  cout << "full integration no log  " << myBB.integratedFlux() << endl;
+  cout << "full integration log     " << myBB.logIntegratedFlux() << endl;
+
+  /*
+  cout << "SZSource " << endl;
+  SzSource mySZ(1);
+  double logSZ = mySZ.logIntegratedFlux(freqmin,freqmax);
+  double SZ    = mySZ.integratedFlux(freqmin,freqmax);
+  cout << "log   " << logSZ << endl;
+  cout << "nolog " << SZ    << endl;
+  cout << "min et max " << mySZ.minFreq() << " :: " << mySZ.maxFreq() << endl;
+  cout << "full integration no log  " << mySZ.integratedFlux() << endl;
+  cout << "full integration log     " << mySZ.logIntegratedFlux() << endl;
+  */
+  cout << "squareFilter " << endl;
+  SquareFilter mySF(0,200);
+  double logSF = mySF.logIntegratedSpect(freqmin,freqmax);
+  double SF    = mySF.IntegratedSpect(freqmin,freqmax);
+  cout << "log   " << logSF << endl;
+  cout << "nolog " << SF    << endl;
+  cout << "min et max " << mySF.minFreq() << " :: " << mySF.maxFreq() << endl;
+  cout << "full integration no log  " << mySF.IntegratedSpect() << endl;
+  cout << "full integration log     " << mySF.logIntegratedSpect() << endl;
+
+  cout << "triangleFilter" << endl;
+  TriangleFilter myTF(10,200,20,40);
+  double logTF = myTF.logIntegratedSpect(10.,freqmax);
+  double TF    = myTF.IntegratedSpect(10.,freqmax);
+  cout << "log   " << logTF << endl;
+  cout << "nolog " << TF    << endl;
+  cout << "min et max " << myTF.minFreq() << " :: " << myTF.maxFreq() << endl;
+  cout << "full integration no log  " << myTF.IntegratedSpect() << endl;
+  cout << "full integration log     " << myTF.logIntegratedSpect() << endl;
+  
+  cout << "....Testing RadSpectraVec..." << endl;
+  int maxVal = 10;
+  Vector vecOfNu(maxVal);
+  Vector vecOfFDeNu(maxVal);
+  for (int i=0; i<maxVal; i++)
+    {
+      vecOfNu(i) = (double)i+30.;
+      vecOfFDeNu(i) = 1.; //*(double)i+1;
+      cout << "entry values  " <<  vecOfNu(i) << " " << vecOfFDeNu(i) << endl;
+    }
+  freqmin = 10.;
+  freqmax = 19.;
+  RadSpectraVec myRSV(vecOfNu, vecOfFDeNu);
+  double RSV    = myRSV.integratedFlux(freqmin,freqmax);
+ 
+  double logRSV = myRSV.logIntegratedFlux(freqmin,freqmax);
+  cout << "nolog int " << RSV    << endl; 
+  cout << "log int   " << logRSV << endl;
+  cout << "min et max " << myRSV.minFreq() << " :: " << myRSV.maxFreq() << endl;
+  cout << "full integration no log  " << myRSV.integratedFlux() << endl;
+  cout << "full integration log     " << myRSV.logIntegratedFlux() << endl;
+
+  cout << "....Testing SpecRespVec..." << endl;
+  int SRVmaxVal = 100;
+  Vector vecOfNuRSV(SRVmaxVal);
+  Vector vecOfFDeNuRSV(SRVmaxVal);
+  for (int i=0; i<SRVmaxVal; i++)
+    {
+      vecOfNuRSV(i) = (double)i+30.;
+      vecOfFDeNuRSV(i) = 10*(double)i+1;
+      //      cout << "entry values  " <<  vecOfNuRSV(i) << " " << vecOfFDeNuRSV(i) << endl;
+    }
+  freqmin = vecOfNuRSV(0);
+  freqmax = vecOfNuRSV(SRVmaxVal-1);
+  SpecRespVec mySRV(vecOfNuRSV, vecOfFDeNuRSV);
+  double SRV    = mySRV.IntegratedSpect(freqmin,freqmax);
+  double logSRV = mySRV.logIntegratedSpect(freqmin,freqmax);
+  cout << "nolog " << SRV    << endl; 
+  cout << "log   " << logSRV << endl;
+  cout << "min et max " << mySRV.minFreq() << " :: " << mySRV.maxFreq() << endl;
+  cout << "full integration no log " << mySRV.IntegratedSpect() << endl;
+  cout << "full integration log    " << mySRV.logIntegratedSpect() << endl;
+  
+  cout << "----------------------------------------------------" << endl;
+  cout << "...Testing filter and integrate !" << endl;
+
+  cout << " BBSF nolog " << myBB.filteredIntegratedFlux(mySF,30.,100.) << endl;
+  cout << " BBSF   log " << myBB.filteredLogIntFlux(mySF,30.,100.) << endl;
+  cout << "full integration no log  " << myBB.integratedFlux(30.,100.) << endl;
+  cout << "full integration log     " << myBB.logIntegratedFlux(30.,100.) << endl;
+  cout << "SF" << endl;
+  cout << "full integration no log  " << mySF.IntegratedSpect(30.,100.) << endl;
+  cout << "full integration log     " << mySF.logIntegratedSpect(30.,100.) << endl;
+
+  cout << " BBTF nolog " << myBB.filteredIntegratedFlux(myTF,30.,100.) << endl;
+  cout << " BBTF   log " << myBB.filteredLogIntFlux(myTF,30.,100.) << endl;
+  cout << "TF" << endl;
+  cout << "full integration no log  " << myTF.IntegratedSpect(30.,100.) << endl;
+  cout << "full integration log     " << myTF.logIntegratedSpect(30.,100.) << endl;
+
+  /*
+  cout << " SZSF nolog " << mySZ.filteredIntegratedFlux(mySF,30.,100.) << endl;
+  cout << " SZSF  log " << mySZ.filteredLogIntFlux(mySF,30.,100.) << endl;
+  cout << "full integration no log  " << mySZ.integratedFlux(30.,100.) << endl;
+  cout << "full integration log     " << mySZ.logIntegratedFlux(30.,100.) << endl;
+
+  cout << " SZTF nolog " << mySZ.filteredIntegratedFlux(myTF,30.,100.) << endl;
+  cout << " SZTF   log " << mySZ.filteredLogIntFlux(myTF,30.,100.) << endl;
+			   */
+  cout << "********" << endl;
+  cout << " RSVSF nolog " << myRSV.filteredIntegratedFlux(mySF,30.,100.) << endl;
+  cout << " RSVSF   log " << myRSV.filteredLogIntFlux(mySF,30.,100.) << endl;
+  cout << "full integration no log  " << myRSV.integratedFlux(30.,100.) << endl;
+  cout << "full integration log     " << myRSV.logIntegratedFlux(30.,100.) << endl;
+
+  cout << " RSVTF nolog " << myRSV.filteredIntegratedFlux(myTF,30.,100.) << endl;
+  cout << " RSVTF   log " << myRSV.filteredLogIntFlux(myTF,30.,100.) << endl;
+  
+
+  return 1;
+} 
Index: /trunk/SophyaLib/SkyT/gaussfilt.cc
===================================================================
--- /trunk/SophyaLib/SkyT/gaussfilt.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/gaussfilt.cc	(revision 601)
@@ -0,0 +1,77 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: gaussfilt.cc,v 1.1.1.1 1999-11-19 16:34:33 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+//---------------
+// C++ Headers --
+//---------------
+#include "machdefs.h"
+#include <iostream.h>
+#include <math.h>
+
+#include "gaussfilt.h"
+
+//----------------
+// Constructor --
+//----------------
+GaussianFilter::GaussianFilter(double nu0, double s, double a, double numin, double numax)
+  : SpectralResponse(numin, numax)
+{
+  if (s < 1.e-19) s = 1.e-19;
+  _s = s;
+  _nu0 = nu0;
+  _a = a;
+}
+
+
+//--------------
+// Destructor --
+//--------------
+GaussianFilter::~GaussianFilter()
+{
+}
+
+//              ---------------------------
+//              --  Function Definitions --
+//              ---------------------------
+
+
+double 
+GaussianFilter::transmission(double nu) const
+{
+  if ((nu < _numin) || (nu > _numax)) return(0.);
+  else { 
+    double tmp = (nu-_nu0)/_s;
+    return(_a * exp(tmp*tmp));
+  }
+}
+
+double 
+GaussianFilter::peakFreq()  const
+{
+return(_nu0);
+}
+
+double 
+GaussianFilter::peakTransmission() const
+{
+return(_a);
+}
+
+void 
+GaussianFilter::Print(ostream& os)  const
+{ 
+  os << "GaussianFilter::Print - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << " T = A * Exp(((nu-nu0)/s)^2) : " << " nu0= " << _nu0 << " sig= " 
+     << _s << "  A= " << _a << endl;  
+  os << "PeakFreq= " << peakFreq() << "  Transmission= " << transmission(peakFreq()) << endl;
+
+
+}
Index: /trunk/SophyaLib/SkyT/gaussfilt.h
===================================================================
--- /trunk/SophyaLib/SkyT/gaussfilt.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/gaussfilt.h	(revision 601)
@@ -0,0 +1,56 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: gaussfilt.h,v 1.1.1.1 1999-11-19 16:34:33 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+#ifndef GAUSSFILT_H_SEEN
+#define GAUSSFILT_H_SEEN
+
+//----------------------
+// Base Class Headers --
+//----------------------
+#include "specresp.h"
+
+//------------------------------------
+// Collaborating Class Declarations --
+//------------------------------------
+
+//              ---------------------
+//              -- Class Interface --
+//              ---------------------
+
+// Spectral response in the form  A Exp(((nu-nu0)/s)^2) 
+
+class GaussianFilter:public SpectralResponse
+{
+public:
+  //Constructors
+  GaussianFilter (double nu0, double s, double a=1., double numin=0., double numax= 9.e49); 
+  
+
+  // destructor
+  virtual ~GaussianFilter();
+
+  //
+  //  Function members
+  //
+
+  virtual double transmission(double) const;
+  virtual double peakFreq() const;
+  virtual double peakTransmission() const;
+
+  virtual void   Print(ostream& os)  const;
+
+protected:
+  double _a, _s, _nu0;
+};
+
+
+#endif
Index: /trunk/SophyaLib/SkyT/nupower.cc
===================================================================
--- /trunk/SophyaLib/SkyT/nupower.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/nupower.cc	(revision 601)
@@ -0,0 +1,55 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: nupower.cc,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//      Aim of the class: To give the energy density for a nupower
+//                        The unity used here is W/m^2/Hz/sr
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+//---------------
+// C++ Headers --
+//---------------
+#include "machdefs.h"
+#include <iostream.h>
+#include <math.h>
+
+#include "nupower.h"
+
+//----------------
+// Constructor --
+//----------------
+PowerLawSpectra::PowerLawSpectra(double a, double b, double nu0, double numin, double numax)
+  : RadSpectra(numin, numax)
+{
+  _a = a;
+  _b = b;
+  _nu0 = nu0;
+}
+
+
+PowerLawSpectra::~PowerLawSpectra()
+{
+}
+
+double
+PowerLawSpectra::flux(double nu) const
+{
+  if ((nu < _numin) || (nu > _numax)) return(0.);
+  else return(  _a * pow((nu-_nu0), _b) );
+}
+
+void
+PowerLawSpectra::Print(ostream& os) const
+{
+  os << "PowerLawSpectra::Print f = a (nu-nu0)^b " << _a << "(nu-" << _nu0 
+     << ") ^ " << _b << endl;
+  os << " - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "MeanFreq= " << meanFreq() << "  Emission= " << flux(meanFreq()) << endl;
+  os << "PeakFreq= " << peakFreq() << "  Emission= " << flux(peakFreq()) << endl;
+
+}
Index: /trunk/SophyaLib/SkyT/nupower.h
===================================================================
--- /trunk/SophyaLib/SkyT/nupower.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/nupower.h	(revision 601)
@@ -0,0 +1,37 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: nupower.h,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+#ifndef NUPOWER_H_SEEN
+#define NUPOWER_H_SEEN
+
+#include "radspec.h"
+#include "convtools.h"
+
+// Power law spectra  f = a (nu-nu0)^b   
+class PowerLawSpectra : public RadSpectra
+{
+public:   //Constructor
+  
+  PowerLawSpectra(double a, double b, double nu0, double numin=0., double numax=9.e49);
+
+  virtual ~PowerLawSpectra();
+
+  virtual double flux(double nu) const;
+
+  // Redefintion of the Print function
+  virtual void      Print(ostream& os)  const;
+
+protected:
+  double _a, _b, _nu0;
+};
+
+#endif
Index: /trunk/SophyaLib/SkyT/radspec.cc
===================================================================
--- /trunk/SophyaLib/SkyT/radspec.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/radspec.cc	(revision 601)
@@ -0,0 +1,236 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: radspec.cc,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//      Aim of the class: To give the energy density
+//                        The unity used here is W/m^2/Hz/sr
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+//---------------
+// C++ Headers --
+//---------------
+#include "machdefs.h"
+#include <iostream.h>
+#include <typeinfo>
+#include <math.h>
+
+#include "radspec.h"
+#include "integ.h"
+
+//----------------
+// Constructor --
+//----------------
+RadSpectra::RadSpectra(double numin, double numax)
+{
+  _numin = numin;
+  _numax = numax;
+}
+
+
+//--------------
+// Destructor --
+//--------------
+RadSpectra::~RadSpectra()
+{
+}
+
+//              ---------------------------
+//              --  Function Definitions --
+//              ---------------------------
+
+double 
+RadSpectra::minFreq()  const
+{
+  return _numin;
+}
+
+double 
+RadSpectra::maxFreq()  const
+{
+  return _numax;
+}
+
+double 
+RadSpectra::meanFreq()  const
+{
+  double result = (_numax+_numin)/2.;
+  return result;
+}
+
+
+
+// peakFreq returns the value of the frequency for the
+// peak of the spectrum. 
+double 
+RadSpectra::peakFreq()  const
+{
+  double maxAnswer = -1.e99;
+  double maxNu = -10;
+  double nu;
+  for (int i=1; i<1000;i++)
+    {
+      nu=(_numax-_numin)*i/1000.+_numin;
+      double lookForMax = flux(nu);
+      if(maxAnswer <= lookForMax) {
+	maxAnswer= lookForMax;
+	maxNu    = nu;
+      }
+    }
+  return maxNu;  
+}
+
+// To change min-max frequency
+void
+RadSpectra::setMinMaxFreq(double numin, double numax)
+{
+  _numin = numin;
+  _numax = numax;
+}
+
+// the RadSpectra_fluxFunction function is used to call TrpzInteg double(double) 
+// (integration over a range of frequencies)
+static RadSpectra* _raypourfinteg = NULL;
+static double RadSpectra_fluxFunction(double nu)
+{
+   return(_raypourfinteg->flux(nu));   
+}
+
+double 
+RadSpectra::integratedFlux(double f1, double f2)  const
+{
+    _raypourfinteg = const_cast<RadSpectra *>(this);
+    TrpzInteg I(RadSpectra_fluxFunction , f1, f2);  
+    _raypourfinteg = NULL;
+    return((double)I);        
+}
+double 
+RadSpectra::integratedFlux()  const
+{
+  return integratedFlux(_numin, _numax);
+}
+
+// integration using the logarithm !!
+// Carefull!! Base 10....
+static RadSpectra* _rayIntLog = NULL;
+
+static double RadSpectra_logFluxFunction(double tau)
+{
+  double value = _rayIntLog->flux(pow(10,tau))*pow(10,tau);
+  return(value);   
+}
+
+double 
+RadSpectra::logIntegratedFlux(double f1, double f2)  const
+{
+  double f1Log = log10(f1);
+  double f2Log = log10(f2);
+  if(f1Log < -1.e99) f1Log = -1.e99;
+  if(f2Log > 1.e99)  f2Log = 1.e99;
+  _rayIntLog = const_cast<RadSpectra *>(this);
+  TrpzInteg I(RadSpectra_logFluxFunction,f1Log,f2Log);
+  double value = (double)I * log(10.);
+  _rayIntLog = NULL;
+  return(value);        
+}
+
+double 
+RadSpectra::logIntegratedFlux()  const
+{
+  return logIntegratedFlux(_numin,_numax);
+}
+
+// the RadSpectra_filteredFlux function is used to call TrpzInteg double(double) 
+// (integration over a range of frequencies with a filter)
+static SpectralResponse* _filter = NULL ;
+static double RadSpectra_filteredFlux(double nu) 
+{
+  double flux = _raypourfinteg->flux(nu);
+  return(flux * _filter->transmission(nu));
+}
+
+
+double 
+RadSpectra::filteredIntegratedFlux(SpectralResponse const& filter, double f1, double f2) const
+{
+  _raypourfinteg = const_cast<RadSpectra *>(this);
+  _filter = const_cast<SpectralResponse *>(&filter);
+  TrpzInteg I(RadSpectra_filteredFlux,f1,f2);  
+  _raypourfinteg = NULL;
+  _filter = NULL;
+  return((double)I);        
+}
+
+double 
+RadSpectra::filteredIntegratedFlux(SpectralResponse const& filter)
+{
+  double minOfMin = filter.minFreq();
+  double maxOfMax = filter.maxFreq();
+  if(minOfMin < this->minFreq()) minOfMin = this->minFreq();
+  if(maxOfMax > this->maxFreq()) maxOfMax = this->maxFreq();  
+  return(filteredIntegratedFlux(filter, minOfMin, maxOfMax ) );
+}
+
+
+// the RadSpectraVec_filteredFlux function is used to call TrpzInteg double(double) 
+// (integration over a range of frequencies with a filter)
+static double RadSpectra_logFilteredFlux(double tau) 
+{
+  double nu = pow(10,tau);
+  double flux = _raypourfinteg->flux(nu)*nu;
+  return(flux * _filter->transmission(nu));
+}
+
+
+double 
+RadSpectra::filteredLogIntFlux(SpectralResponse const& filter, double f1, double f2) const
+{
+
+//   double f1Log = log10(f1);
+//   double f2Log = log10(f2);
+//   if(f1Log < -1.e99) f1Log = -1.e99;
+//   if(f2Log > 1.e99)  f2Log = 1.e99;
+//   _rayIntLog = const_cast<RadSpectra *>(this);
+//   TrpzInteg I(RadSpectra_logFluxFunction,f1Log,f2Log);
+//   double value = (double)I * log(10.);
+//   return(value);        
+//   _rayIntLog = NULL;
+  
+   _raypourfinteg = NULL;
+   _filter = NULL;
+   double f1Log = log10(f1);
+   double f2Log = log10(f2);
+   if(f1Log < -1.e99) f1Log = -1.e99;
+   if(f2Log > 1.e99)  f2Log = 1.e99;
+   _raypourfinteg = const_cast<RadSpectra *>(this);
+   _filter = const_cast<SpectralResponse *>(&filter);
+   TrpzInteg I(RadSpectra_logFilteredFlux,f1Log,f2Log);  
+   _raypourfinteg = NULL;
+   _filter = NULL;
+   return((double)I* log(10.));        
+}
+
+double 
+RadSpectra::filteredLogIntFlux(SpectralResponse const& filter)
+{
+  return(filteredLogIntFlux(filter, filter.minFreq(), filter.maxFreq() ) );
+}
+
+
+
+void
+RadSpectra::Print(ostream& os) const
+{
+  //  os << "RadSpectra::Print (" << typeid(*this).name() 
+  //     << ") - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "RadSpectra::Print  - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "MeanFreq= " << meanFreq() << "  Emission= " << flux(meanFreq()) << endl;
+  os << "PeakFreq= " << peakFreq() << "  Emission= " << flux(peakFreq()) << endl;
+
+}
+
+
Index: /trunk/SophyaLib/SkyT/radspec.h
===================================================================
--- /trunk/SophyaLib/SkyT/radspec.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/radspec.h	(revision 601)
@@ -0,0 +1,89 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: radspec.h,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+#ifndef RADSPEC_H_SEEN
+#define RADSPEC_H_SEEN
+
+#include "machdefs.h"
+#include <iostream.h>
+
+#include "specresp.h"
+
+//              ---------------------
+//              -- Class Interface --
+//              ---------------------
+
+// ***IMPORTANT*** All frequencies are expressed in GHz  (10^9 Hz)
+//                 and flux in Watt/m^2/sr
+
+class RadSpectra
+{
+
+public:
+  //Constructor
+  RadSpectra (double numin=0., double numax=9.e49);
+
+  // destructor
+  virtual ~RadSpectra();
+
+  //
+  //  Member Functions
+  //
+
+  //  The flux() function is virtual: 
+  //  It computes the flux per units of W/m^2/Hz/st
+  // implemented in the subclasses
+  virtual double flux(double nu) const =0;
+  inline double operator() (double nu) const { return(flux(nu)); }
+
+  virtual double minFreq() const;
+  virtual double maxFreq() const;
+  virtual double meanFreq() const;
+  virtual double peakFreq() const;
+
+  // For changing the min/max frequency
+  virtual void	 setMinMaxFreq(double numin, double numax);
+
+  // The IntegratedFlux() function
+  // integrates the flux function between two frequencies: 
+  // numin and numax (in GHz)
+  virtual double integratedFlux( double numin, double numax) const ;
+  virtual double logIntegratedFlux( double numin, double numax) const ;
+
+  // same but integrated over the full spectrum of frequencies
+  virtual  double integratedFlux() const ;
+  virtual  double logIntegratedFlux() const ;
+  
+  // filteredIntegratedFlux integrates the flux function while
+  // taking into account the spectralResponse of a filter
+  virtual double filteredIntegratedFlux(SpectralResponse const & filter, 
+					double numin, double numax) const ;
+  virtual double filteredLogIntFlux(SpectralResponse const & filter, 
+					double numin, double numax) const ;
+
+  // same but integrates over the filter's frequency band
+  virtual double filteredIntegratedFlux(SpectralResponse const & filter);
+  virtual double filteredLogIntFlux(SpectralResponse const & filter);
+
+  virtual void      Print(ostream& os)  const;
+
+protected:
+  double _numin;
+  double _numax;
+};
+
+// Definition de l'operateur <<
+inline ostream& operator << (ostream& s,  RadSpectra const & rs)
+  {  rs.Print(s);  return(s);  }
+
+
+#endif
Index: /trunk/SophyaLib/SkyT/radspecvector.cc
===================================================================
--- /trunk/SophyaLib/SkyT/radspecvector.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/radspecvector.cc	(revision 601)
@@ -0,0 +1,97 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: radspecvector.cc,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//      Aim of the class: To give the energy density
+//                        The unity used here is W/m^2/Hz/sr
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+//---------------
+// C++ Headers --
+//---------------
+#include "machdefs.h"
+#include <iostream.h>
+#include <math.h>
+// #include <typeinfo>
+
+#include "radspecvector.h"
+#include "pexceptions.h"
+
+//----------------
+// Constructor --
+//----------------
+RadSpectraVec::RadSpectraVec(Vector const & nu, Vector const & fdenu, double numin, double numax)
+        : RadSpectra(numin, numax)
+{
+  if(nu.NElts() != fdenu.NElts())  
+    throw SzMismatchError("RadSpectraVec::RadSpectraVec() - Non equal vector sizes");
+  _vecOfNu = nu;
+  _vecOfFDeNu = fdenu;
+  _numin = nu(0);
+  _numax = nu(nu.NElts()-1);
+  if(_vecOfNu.NElts() != _vecOfFDeNu.NElts()) cout << "vectors are not compatible" << exit;
+}
+
+
+//--------------
+// Destructor --
+//--------------
+RadSpectraVec::~RadSpectraVec()
+{
+}
+
+//              ---------------------------
+//              --  Function Definitions --
+//              ---------------------------
+
+
+
+
+double 
+RadSpectraVec::flux(double nu) const
+{
+  if ( (nu < _numin) || (nu > _numax) ) return(0.);
+  double value = -99;
+  int sizeVecOfNu = _vecOfNu.NElts();
+  double minVal = _vecOfNu(0);
+  if(nu <= minVal) return _vecOfFDeNu(0);
+  if(nu >=  _vecOfNu(sizeVecOfNu-1)) return _vecOfFDeNu(sizeVecOfNu-1);
+  
+  for (int i=1; i<sizeVecOfNu; i++)
+    {
+      if(nu>= minVal && nu < _vecOfNu(i))
+	{
+	  double up = _vecOfFDeNu(i) ;
+	  double down = _vecOfFDeNu(i-1);
+	  double xmin = _vecOfNu(i-1);
+	  double xmax = _vecOfNu(i);
+	  double a = ((up-down)/(xmax-xmin));
+	  value = a*nu+(up-a*xmax);
+	  return value;
+	}
+      else 
+	{
+	  minVal = _vecOfNu(i);
+	}	
+    }
+  return value;
+}
+
+
+void
+RadSpectraVec::Print(ostream& os) const
+{
+  //  os << "RadSpectraVec::Print (" << typeid(*this).name() 
+  //     << ") - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "RadSpectraVec::Print  - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "MeanFreq= " << meanFreq() << "  Emission= " << flux(meanFreq()) << endl;
+  os << "PeakFreq= " << peakFreq() << "  Emission= " << flux(peakFreq()) << endl;
+
+}
+
+
Index: /trunk/SophyaLib/SkyT/radspecvector.h
===================================================================
--- /trunk/SophyaLib/SkyT/radspecvector.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/radspecvector.h	(revision 601)
@@ -0,0 +1,54 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+#ifndef RADSPECTRAVEC_H_SEEN
+#define RADSPECTRAVEC_H_SEEN
+#include "machdefs.h"
+#include <iostream.h>
+
+#include "radspec.h"
+#include "tvector.h"
+
+//              ---------------------
+//              -- Class Interface --
+//              ---------------------
+
+// ***IMPORTANT*** All frequencies are expressed in GHz  (10^9 Hz)
+//                 and flux in Watt/m^2/sr
+
+class RadSpectraVec : public RadSpectra
+{
+
+public:
+  //Constructor
+  RadSpectraVec (Vector const & nu, Vector const & fdenu, double numin=0., double numax=9.e49);
+
+  // destructor
+  virtual ~RadSpectraVec();
+
+  //
+  //  Member Functions
+  //
+
+  //  The flux() function is virtual: 
+  //  It computes the flux per units of W/m^2/Hz/st
+  // implemented in the subclasses
+  virtual double flux(double nu) const;
+
+  virtual void   Print(ostream& os)  const;
+
+protected:
+  Vector _vecOfNu;
+  Vector _vecOfFDeNu;
+  
+};
+
+#endif
Index: /trunk/SophyaLib/SkyT/skm.d
===================================================================
--- /trunk/SophyaLib/SkyT/skm.d	(revision 601)
+++ /trunk/SophyaLib/SkyT/skm.d	(revision 601)
@@ -0,0 +1,13 @@
+#  SKYMIX  NbComponents  M_Gorski
+@SKYMIX  2  64
+#  FILTER  Nu0 Sigma_Nu Tmax NuMin NuMax
+@FILTER 143 12 0.84 50. 300.
+
+#  Sky components 
+#  SCMFILE  FITSfilename [Normalisation]
+@SCMFILE1  comp1.fits  
+@SCMFILE2  comp2.fits    0.25 
+
+#  SCSPEC1  FITSfilename  Fmin Fmax 
+@SCSPEC1  spec1.fits  
+@SCSPEC2  spec2.fits
Index: /trunk/SophyaLib/SkyT/skymixer.cc
===================================================================
--- /trunk/SophyaLib/SkyT/skymixer.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/skymixer.cc	(revision 601)
@@ -0,0 +1,159 @@
+#include "machdefs.h"
+#include <stdlib.h>
+#include <stdio.h>
+
+#include <iostream.h>
+#include <math.h>
+
+#include <string>
+#include <vector>
+
+#include "timing.h"
+#include "sambainit.h"
+#include "pexceptions.h"
+#include "datacards.h"
+#include "fitsioserver.h"
+
+#include "spheregorski.h"
+
+#include "radspecvector.h"
+#include "blackbody.h"
+#include "nupower.h"
+
+#include "squarefilt.h"
+#include "trianglefilt.h"
+#include "specrespvector.h"
+#include "gaussfilt.h"
+
+// ------ Function declaration 
+void addComponent(SpectralResponse&  sr, PixelMap<float>& finalMap, 
+                  PixelMap<float>& mapToAdd, RadSpectra& rs, double K=1.);
+void addComponent(SpectralResponse&  sr, PixelMap<double>& finalMap, 
+                  PixelMap<double>& mapToAdd, RadSpectra& rs, double K=1.);
+
+// -------------------------------------------------------------------------
+//                             main program 
+// -------------------------------------------------------------------------
+int main(int narg, char * arg[])
+{
+  if ((narg < 2) || (strcmp(arg[1], "-h") == 0) ) {
+    cout << " skymixer / Error args \n Usage: skymixer parameterFile " << endl;
+    exit(0);
+  }
+
+InitTim();
+
+string msg;
+int rc = 0;
+  
+try {
+  string dcard = arg[1];
+  cout << " Decoding parameters from file " << dcard << endl;
+  DataCards dc(dcard);
+
+  //  Cheking datacards
+  if ( (!dc.HasKey("SKYMIX")) || (!dc.HasKey("FILTER")) ) {
+     rc = 71; 
+     msg = "Invalid parameters  - NO @SKYMIX or @FILTER card ";
+     goto problem;
+  }
+
+  // Decoding number of component and pixelisation parameter
+  int mg = 32;
+  int ncomp = 0;
+  ncomp = dc.IParam("SKYMIX", 0, 0);
+  mg = dc.IParam("SKYMIX", 1, 32);
+  if (ncomp < 1)  {
+    msg = "Invalid parameters  - Check datacards @SKYMIX ";
+    rc = 72;
+    goto problem;
+  }
+
+  int kc;
+  string key;
+  char buff[256];
+  bool pb = false;
+  for(kc=0; kc<ncomp; kc++) {
+    sprintf(buff, "SCMFILE%d", kc+1);
+    key = buff;
+    if (dc.NbParam(key) < 1) {   
+      msg = "Missing or invalid card : " + key;
+      pb = true;  break;
+    }
+    sprintf(buff, "SCSPEC%d", kc+1);
+    key = buff;
+    if (dc.NbParam(key) < 1)  {   
+      msg = "Missing or invalid card : " + key;
+      pb = true;  break;
+    }
+
+    }
+
+  if (pb)  {
+    rc = 72;
+    goto problem;
+  }
+
+  cout << " skymix/Info : NComp = " <<  ncomp << " M_SphereGorski= " << mg << endl;
+
+  SphereGorski<float> outgs(mg);
+  cout << " Creating Output Gorski Map  NbPixels= " << outgs.NbPixels() << endl;
+  outgs.SetPixels(0.);
+
+  // Decoding detection pass-band filter 
+  double nu0 = dc.DParam("FILTER", 0, 10.);
+  double s = dc.DParam("FILTER", 1, 1.);
+  double a = dc.DParam("FILTER", 2, 1.);
+  double numin = dc.DParam("FILTER", 3, 0.1);
+  double numax = dc.DParam("FILTER", 4, 9999);
+  GaussianFilter filt(nu0, s, a, numin, numax);
+  cout << " Filter decoded - Created " << endl;
+  cout << filt << endl;
+
+  PrtTim(" After FilterCreation ");
+  }
+
+catch (PException exc) {
+  msg = exc.Msg();
+  cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
+  rc = 70;
+  }   
+
+
+problem:
+if (rc == 0)  return(0);
+cerr << " Error condition -> Rc= " << rc << endl;
+cerr << " Msg= " << msg << endl;
+return(rc);
+}
+
+// template <class T> 
+void addComponent(SpectralResponse&  sr, PixelMap<float>& finalMap, 
+                  PixelMap<float>& mapToAdd, RadSpectra& rs, double K)
+{
+  // finalMap = finalMap + coeff* mapToAdd
+  // coeff    =  convolution of sr and rs 
+  // compute the coefficient corresponding to mapToAdd
+  if (finalMap.NbPixels() != mapToAdd.NbPixels())   
+    throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
+  double coeff = rs.filteredIntegratedFlux(sr) * K;
+  for(int i=0; i<finalMap.NbPixels(); i++)
+    {
+      finalMap(i) += coeff * mapToAdd(i);
+    }
+}
+
+void addComponent(SpectralResponse&  sr, PixelMap<double>& finalMap, 
+                  PixelMap<double>& mapToAdd, RadSpectra& rs, double K)
+{
+  // finalMap = finalMap + coeff* mapToAdd
+  // coeff    =  convolution of sr and rs 
+  // compute the coefficient corresponding to mapToAdd
+  if (finalMap.NbPixels() != mapToAdd.NbPixels())   
+    throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
+  double coeff = rs.filteredIntegratedFlux(sr) * K;
+  for(int i=0; i<finalMap.NbPixels(); i++)
+    {
+      finalMap(i) += coeff * mapToAdd(i);
+    }
+}
Index: /trunk/SophyaLib/SkyT/specresp.cc
===================================================================
--- /trunk/SophyaLib/SkyT/specresp.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/specresp.cc	(revision 601)
@@ -0,0 +1,161 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: specresp.cc,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//      Aim of the class: To give the energy density
+//                        The unity used here is W/m^2/Hz/sr
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+//---------------
+// C++ Headers --
+//---------------
+#include <typeinfo>
+#include <math.h>
+#include "specresp.h"
+#include "integ.h"
+
+//----------------
+// Constructor --
+//----------------
+SpectralResponse::SpectralResponse(double numin, double numax)
+{
+  _numin = numin;
+  _numax = numax;
+}
+
+
+//--------------
+// Destructor --
+//--------------
+SpectralResponse::~SpectralResponse()
+{
+}
+
+//              ---------------------------
+//              --  Function Definitions --
+//              ---------------------------
+
+
+double 
+SpectralResponse::minFreq()  const 
+{
+  return _numin;
+}
+
+double 
+SpectralResponse::maxFreq()  const
+{
+  return _numax;
+}
+
+double 
+SpectralResponse::meanFreq()  const
+{
+  return (_numax+_numin)/2.;
+}
+
+ 
+
+// peakFreq returns the value of the frequency for the
+// peak of the spectrum. 
+double 
+SpectralResponse::peakFreq()  const
+{
+  double maxAnswer = -1.e99;
+  double maxNu = -10;
+  double nu;
+  for (int i=1; i<1000;i++)
+    {
+      nu=(_numax-_numin)*i/1000.+_numin;
+      double lookForMax =transmission(nu);
+      if(maxAnswer <= lookForMax) {
+	maxAnswer= lookForMax;
+	maxNu    = nu;
+      }
+    }
+  return maxNu;    
+}
+
+
+double 
+SpectralResponse::peakTransmission()  const
+{
+  double nuPeak = this->peakFreq();
+  return transmission(nuPeak);
+}
+
+static SpectralResponse* _mySpecResp = NULL;
+
+static double SpectralResponse_transmission(double nu)
+{
+   return(_mySpecResp->transmission(nu));   
+}
+
+
+double 
+SpectralResponse::IntegratedSpect()  const
+{
+  double value = this->IntegratedSpect(_numin, _numax);
+  return value;
+}
+
+double 
+SpectralResponse::IntegratedSpect(double numin, double numax) const 
+{
+  _mySpecResp = const_cast<SpectralResponse *>(this);
+  if(numin <= _numin) numin = _numin;
+  if(numax >= _numax) numax = _numax;
+  TrpzInteg I(SpectralResponse_transmission , numin, numax);  
+  return((double)I);
+  _mySpecResp= NULL;
+}
+
+
+static SpectralResponse* _myLogSpecResp = NULL;
+
+static double SpectralResponse_logTransmission(double tau)
+{
+  double value =  _myLogSpecResp->transmission(pow(10,tau))*pow(10,tau);
+  return(value);   
+}
+
+double 
+SpectralResponse::logIntegratedSpect(double numin, double numax) const 
+{
+  if(numin <= _numin) numin = _numin;
+  if(numax >= _numax) numax = _numax;
+  if(numin == 0) numin = 1.e-99;
+  double f1Log = log10(numin);
+  double f2Log = log10(numax);
+  if(f1Log < -1.e99) f1Log = -1.e99;
+  if(f2Log > 1.e99)  f2Log = 1.e99;
+  _myLogSpecResp = const_cast<SpectralResponse *>(this); 
+  TrpzInteg I(SpectralResponse_logTransmission ,f1Log,f2Log); 
+  return((double)I*log(10.));
+  _myLogSpecResp= NULL;
+}
+
+double 
+SpectralResponse::logIntegratedSpect()  const
+{
+  double value = this->logIntegratedSpect(_numin, _numax);
+  return value;
+}
+
+
+void
+SpectralResponse::Print(ostream& os) const 
+{
+
+  //   os << "SpectralResponse::Print (" << typeid(*this).name() 
+  //<< ") - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "SpectralResponse::Print - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "MeanFreq= " << meanFreq() << "  Transmission= " << transmission(meanFreq()) << endl;
+  os << "PeakFreq= " << peakFreq() << "  Transmission= " << transmission(peakFreq()) << endl;
+
+}
Index: /trunk/SophyaLib/SkyT/specresp.h
===================================================================
--- /trunk/SophyaLib/SkyT/specresp.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/specresp.h	(revision 601)
@@ -0,0 +1,70 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: specresp.h,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+#ifndef SPECRESP_H_SEEN
+#define SPECRESP_H_SEEN
+
+#include "machdefs.h"
+#include <iostream.h>
+
+//------------------------------------
+// Collaborating Class Declarations --
+//------------------------------------
+
+//              ---------------------
+//              -- Class Interface --
+//              ---------------------
+
+// ***IMPORTANT*** All frequencies are expressed in GHz  (10^9 Hz)
+
+class SpectralResponse
+{
+
+public:
+  //Constructor
+  SpectralResponse(double numin=0., double numax=1.);
+
+  // destructor
+  virtual ~SpectralResponse();
+
+  //
+  //  Member Functions
+  //
+
+  //  The transmission() function is virtual: 
+  virtual double transmission(double nu) const =0 ; 
+
+  inline double operator() (double nu) const { return(transmission(nu)); }
+
+  virtual void      Print(ostream& os)  const;
+
+  virtual double minFreq() const;
+  virtual double maxFreq() const;
+  virtual double meanFreq() const;
+  virtual double peakFreq() const;
+  virtual double peakTransmission() const;
+  virtual double IntegratedSpect(double numin, double numax) const ;
+  virtual double logIntegratedSpect(double numin, double numax) const ;
+  virtual double IntegratedSpect() const ;
+  virtual double logIntegratedSpect() const ;
+
+protected:
+  double _numin;
+  double _numax;
+};
+
+
+// definition of the << operator
+inline ostream& operator << (ostream& s,  SpectralResponse const & sr)
+  {  sr.Print(s);  return(s);  }
+
+
+#endif
Index: /trunk/SophyaLib/SkyT/specrespvector.cc
===================================================================
--- /trunk/SophyaLib/SkyT/specrespvector.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/specrespvector.cc	(revision 601)
@@ -0,0 +1,97 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: specrespvector.cc,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//      Aim of the class: To give the energy density
+//                        The unity used here is W/m^2/Hz/sr
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+//---------------
+// C++ Headers --
+//---------------
+#include "machdefs.h"
+#include <iostream.h>
+#include <math.h>
+// #include <typeinfo>
+
+#include "specrespvector.h"
+#include "pexceptions.h"
+
+//----------------
+// Constructor --
+//----------------
+SpecRespVec::SpecRespVec(Vector const & nu, Vector const & fdenu, double numin, double numax)
+        : SpectralResponse(numin, numax)
+{
+  if(nu.NElts() != fdenu.NElts())  
+    throw SzMismatchError("SpecRespVec::SpecRespVec() - Non equal vector sizes");
+  _vecOfNu = nu;
+  _vecOfFDeNu = fdenu;
+  _numin = nu(0);
+  _numax = nu(nu.NElts()-1);
+  if(_vecOfNu.NElts() != _vecOfFDeNu.NElts()) cout << "vectors are not compatible" << exit;
+}
+
+
+//--------------
+// Destructor --
+//--------------
+SpecRespVec::~SpecRespVec()
+{
+}
+
+//              ---------------------------
+//              --  Function Definitions --
+//              ---------------------------
+
+
+
+ 
+double 
+SpecRespVec::transmission(double nu)  const 
+{
+  if ( (nu < _numin) || (nu > _numax) ) return(0.);
+  double value = -99;
+  int sizeVecOfNu = _vecOfNu.NElts();
+  double minVal = _vecOfNu(0);
+  if(nu <= minVal) return _vecOfFDeNu(0);
+  if(nu >=  _vecOfNu(sizeVecOfNu-1)) return _vecOfFDeNu(sizeVecOfNu-1);
+
+  for (int i=1; i<sizeVecOfNu; i++)
+    {
+      if(nu>= minVal && nu < _vecOfNu(i))
+	{
+	  double up = _vecOfFDeNu(i) ;
+	  double down = _vecOfFDeNu(i-1);
+	  double xmin = _vecOfNu(i-1);
+	  double xmax = _vecOfNu(i);
+	  double a = ((up-down)/(xmax-xmin));
+	  value = a*nu+(up-a*xmax);
+	  return value;
+	}
+      else 
+	{
+	  minVal = _vecOfNu(i);
+	}	
+    }
+  return value;
+}
+
+
+
+void
+SpecRespVec::Print(ostream& os) const 
+{
+
+  //   os << "SpecRespVec::Print (" << typeid(*this).name() 
+  //<< ") - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "SpecRespVec ::Print - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
+  os << "MeanFreq= " << meanFreq() << "  Transmission= " << transmission(meanFreq()) << endl;
+  os << "PeakFreq= " << peakFreq() << "  Transmission= " << transmission(peakFreq()) << endl;
+
+}
Index: /trunk/SophyaLib/SkyT/specrespvector.h
===================================================================
--- /trunk/SophyaLib/SkyT/specrespvector.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/specrespvector.h	(revision 601)
@@ -0,0 +1,57 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: specrespvector.h,v 1.1.1.1 1999-11-19 16:34:32 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+#ifndef SPECRESPVEC_H_SEEN
+#define SPECRESPVEC_H_SEEN
+
+#include "machdefs.h"
+#include <iostream.h>
+
+#include "specresp.h"
+#include "tvector.h"
+
+//------------------------------------
+// Collaborating Class Declarations --
+//------------------------------------
+
+//              ---------------------
+//              -- Class Interface --
+//              ---------------------
+
+// ***IMPORTANT*** All frequencies are expressed in GHz  (10^9 Hz)
+
+class SpecRespVec :  public SpectralResponse
+{
+
+public:
+  //Constructor
+  SpecRespVec(Vector const & nu, Vector const & fdenu, double numin=0., double numax=1.);
+
+  // destructor
+  virtual ~SpecRespVec();
+
+  //
+  //  Member Functions
+  //
+
+  //  The transmission() function is virtual: 
+  virtual double transmission(double nu) const; 
+
+  virtual void      Print(ostream& os)  const;
+
+
+protected:
+  Vector _vecOfNu;
+  Vector _vecOfFDeNu;
+};
+
+
+#endif
Index: /trunk/SophyaLib/SkyT/squarefilt.cc
===================================================================
--- /trunk/SophyaLib/SkyT/squarefilt.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/squarefilt.cc	(revision 601)
@@ -0,0 +1,63 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: squarefilt.cc,v 1.1.1.1 1999-11-19 16:34:33 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+//---------------
+// C++ Headers --
+//---------------
+#include "machdefs.h"
+#include <iostream.h>
+
+#include "squarefilt.h"
+
+//----------------
+// Constructor --
+//----------------
+SquareFilter::SquareFilter(double numin, double numax)
+  : SpectralResponse(numin, numax)
+{
+  _nuPeak = (numin+numax)/2.;
+  _peakTransmission = transmission(_nuPeak);
+}
+
+
+//--------------
+// Destructor --
+//--------------
+SquareFilter::~SquareFilter()
+{
+}
+
+//              ---------------------------
+//              --  Function Definitions --
+//              ---------------------------
+
+
+double 
+SquareFilter::transmission(double nu) const
+{
+  if(nu < -1.e99) nu = -1.e99;
+  if(nu > 1.e99) nu = 1.e99;
+  
+  if(nu>=_numin && nu<=_numax) return 1.;
+  return 0.;
+}
+
+double 
+SquareFilter::peakFreq() const 
+{
+return _nuPeak;
+}
+
+double 
+SquareFilter::peakTransmission() const 
+{
+return _peakTransmission;
+}
Index: /trunk/SophyaLib/SkyT/squarefilt.h
===================================================================
--- /trunk/SophyaLib/SkyT/squarefilt.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/squarefilt.h	(revision 601)
@@ -0,0 +1,53 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: squarefilt.h,v 1.1.1.1 1999-11-19 16:34:33 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+#ifndef SQUAREFIL_H_SEEN
+#define SQUAREFIL_H_SEEN
+
+//----------------------
+// Base Class Headers --
+//----------------------
+#include "specresp.h"
+
+//------------------------------------
+// Collaborating Class Declarations --
+//------------------------------------
+
+//              ---------------------
+//              -- Class Interface --
+//              ---------------------
+
+class SquareFilter:public SpectralResponse
+{
+public:
+  //Constructors
+  SquareFilter (double numin, double numax); 
+  
+
+  // destructor
+  virtual ~SquareFilter();
+
+  //
+  //  Function members
+  //
+
+  virtual double transmission(double) const;
+  virtual double peakFreq() const;
+  virtual double peakTransmission() const;
+  
+protected:
+  double _peakTransmission;
+  double _nuPeak;
+};
+
+
+#endif
Index: /trunk/SophyaLib/SkyT/trianglefilt.cc
===================================================================
--- /trunk/SophyaLib/SkyT/trianglefilt.cc	(revision 601)
+++ /trunk/SophyaLib/SkyT/trianglefilt.cc	(revision 601)
@@ -0,0 +1,62 @@
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: trianglefilt.cc,v 1.1.1.1 1999-11-19 16:34:33 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+
+//---------------
+// C++ Headers --
+//---------------
+#include "machdefs.h"
+#include <iostream.h>
+
+#include "trianglefilt.h"
+
+//----------------
+// Constructor --
+//----------------
+TriangleFilter::TriangleFilter(double numin, double numax,double nupeak, double valueatpeak)
+   : SpectralResponse(numin, numax)
+{
+  _nuPeak = nupeak ;
+  _peakTransmission = valueatpeak;
+}
+
+//--------------
+// Destructor --
+//--------------
+TriangleFilter::~TriangleFilter()
+{
+}
+
+//              ---------------------------
+//              --  Function Definitions --
+//              ---------------------------
+
+double 
+TriangleFilter::transmission(double nu) const
+{
+  if(nu < -1.e99) nu = -1.e99;
+  if(nu > 1.e99) nu = 1.e99;
+  
+  if(nu>_numin && nu<_nuPeak) return (_peakTransmission*nu/_nuPeak);
+  if(nu>=_nuPeak && nu<_numax) return (_peakTransmission*nu/(_nuPeak-_numax)+_peakTransmission-_peakTransmission*_nuPeak/(_nuPeak-_numax));
+  return 0.;
+}
+
+double 
+TriangleFilter::peakFreq() const 
+{
+return _nuPeak;
+}
+
+double 
+TriangleFilter::peakTransmission() const 
+{
+return _peakTransmission;
+}
Index: /trunk/SophyaLib/SkyT/trianglefilt.h
===================================================================
--- /trunk/SophyaLib/SkyT/trianglefilt.h	(revision 601)
+++ /trunk/SophyaLib/SkyT/trianglefilt.h	(revision 601)
@@ -0,0 +1,53 @@
+// This may look like C code, but it is really -*- C++ -*-
+//--------------------------------------------------------------------------
+// File and Version Information:
+//      $Id: trianglefilt.h,v 1.1.1.1 1999-11-19 16:34:33 ansari Exp $
+//
+// Description:
+//
+// History (add to end):
+//      Sophie   Oct, 1999  - creation
+//
+//------------------------------------------------------------------------
+#ifndef TRIANGLEFIL_H_SEEN
+#define TRIANGLEFIL_H_SEEN
+
+//----------------------
+// Base Class Headers --
+//----------------------
+#include "specresp.h"
+
+//------------------------------------
+// Collaborating Class Declarations --
+//------------------------------------
+
+//              ---------------------
+//              -- Class Interface --
+//              ---------------------
+
+class TriangleFilter:public SpectralResponse
+{
+public:
+  //Constructors
+  TriangleFilter (double numin, double numax,double nupeak, double valueatpeak); 
+  
+
+  // destructor
+  ~TriangleFilter();
+
+  //
+  //  Function members
+  //
+
+  virtual double transmission(double) const;  
+  virtual double peakFreq() const; 
+  virtual double peakTransmission() const;
+
+protected:
+  double _peakTransmission;
+  double _nuPeak;
+};
+
+
+
+#endif
