double uflux(flux,wave)
double *flux, *wave;
{
  static double uwave[] = {3150.0,3175.0,3200.0,3225.0,3250.0,3275.0,3300.0,
	    3325.0,3350.0,3375.0,3400.0,3425.0,3450.0,3475.0,3500.0,3525.0,
	    3550.0,3575.0,3600.0,3625.0,3650.0,3675.0,3700.0,3725.0,3750.0,
	    3775.0,3800.0,3825.0};
  static double utrans[] = {0.00,0.60,5.92,11.80,17.80,23.80,29.12,33.32,
       36.08,38.12,39.28,39.44,39.08,38.52,37.08,35.32,33.08,30.08,26.16,
       21.68,16.92,11.76,7.68,4.68,2.16,0.88,0.40,0.00};
  double Wave,y1,y2,ut,integral=0.0,a1,a2,A=0.0;
  int i;
  int k = 1;
  int j=0;

  y1 = 0.0;
  a1 = 0.0;
  j = 0;
  while(wave[j] < uwave[27]) {
    if(wave[j] <= uwave[0]) {
      j++;
      continue;
    }
    if(wave[j] >= uwave[27]) break;
    for(i=k;i<=27;i++) {
      if(wave[j] > uwave[i-1] && wave[j] <= uwave[i]) {
	k = i;
	break;
      }
    }
    ut =  utrans[k-1] + (utrans[k]-utrans[k-1])*(wave[j]-uwave[k-1])/
	  (uwave[k]-uwave[k-1]);
    y2 = ut*flux[j]*ip21(wave[j])*bk7(wave[j])*aluminum(wave[j])*
	 aluminum(wave[j])*atm(wave[j],1.3);
    a2 = ut*ip21(wave[j])*bk7(wave[j])*aluminum(wave[j])*aluminum(wave[j])
         *atm(wave[j],1.3);
    integral += 0.5*(y1+y2)*(wave[j]-wave[j-1]);
    A += 0.5*(a1+a2)*(wave[j]-wave[j-1]);
    a1 = a2;
    y1 = y2;
    j++;
  }
  return(integral/A);
}

double vflux(flux,wave)
double *flux,*wave;
{
  static double vwave[] = {3750.0,3775.0,3800.0,3825.0,3850.0,3875.0,3900.0,
     3925.0,3950.0,3975.0,4000.0,4025.0,4050.0,4075.0,4100.0,4125.0,4150.0,
     4175.0,4200.0,4225.0,4250.0,4275.0,4300.0,4325.0,4350.0,4375.0,4400.0,
     4425.0,4450.0};
  static double vtrans[] = {0.00,0.16,0.36,0.84,1.48,2.28,3.04,4.80,7.84,
	13.00,20.00,29.88,40.00,47.28,49.32,48.00,43.52,37.20,28.20,18.12,
	11.12,6.72,4.00,2.68,2.00,1.40,0.72,0.36,0.00};
  double Wave,y1,y2,vt,integral=0.0,A=0.0,a1,a2;
  int i;
  int k = 1;
  int j=0;

  y1 = 0.0;
  a1 = 0.0;
  j = 0;
  while(wave[j] < vwave[28]) {
    if(wave[j] <= vwave[0]) {
      j++;
      continue;
    }
    if(wave[j] >= vwave[28]) break;
    for(i=k;i<=28;i++) {
      if(wave[j] > vwave[i-1] && wave[j] <= vwave[i]) {
	k = i;
	break;
      }
    }
    vt =  vtrans[k-1] + (vtrans[k]-vtrans[k-1])*(wave[j]-vwave[k-1])/
	  (vwave[k]-vwave[k-1]);
    y2 = vt*flux[j]*ip21(wave[j])*bk7(wave[j])*aluminum(wave[j])*
	 aluminum(wave[j])*atm(wave[j],1.3);
    a2 = vt*ip21(wave[j])*bk7(wave[j])*aluminum(wave[j])*aluminum(wave[j])*
         atm(wave[j],1.3);
    integral += 0.5*(y1+y2)*(wave[j]-wave[j-1]);
    A += 0.5*(a1+a2)*(wave[j]-wave[j-1]);
    a1 = a2;
    y1 = y2;
    j++;
  }
  return(integral/A);
}


double bflux(flux,wave)
double *flux,*wave;
{
  static double bwave[] = {4350.0,4375.0,4400.0,4425.0,4450.0,4475.0,4500.0,
     4525.0,4550.0,4575.0,4600.0,4625.0,4650.0,4675.0,4700.0,4725.0,4750.0,
     4775.0,4800.0,4825.0,4850.0,4875.0,4900.0,4925.0,4950.0,4975.0,5000.0,
     5025.0,5050.0};
  static double btrans[] = {0.00,0.40,0.96,1.60,2.32,3.60,5.00,8.00,12.40,
    20.00,30.20,40.40,45.88,46.88,45.04,38.20,27.84,17.40,11.00,6.88,3.92,
    2.52,1.52,1.24,0.92,0.72,0.40,0.20,0.00};
  double Wave,y1,y2,bt,integral=0.0,A=0.0,a1,a2;
  int i;
  int k = 1;
  int j=0;

  y1 = 0.0;
  a1 = 0.0;
  j = 0;
  while(wave[j] < bwave[28]) {
    if(wave[j] <= bwave[0]) {
      j++;
      continue;
    }
    if(wave[j] >= bwave[28]) break;
    for(i=k;i<=28;i++) {
      if(wave[j] > bwave[i-1] && wave[j] <= bwave[i]) {
	k = i;
	break;
      }
    }
    bt =  btrans[k-1] + (btrans[k]-btrans[k-1])*(wave[j]-bwave[k-1])/
	  (bwave[k]-bwave[k-1]);
    y2 = bt*flux[j]*ip21(wave[j])*bk7(wave[j])*aluminum(wave[j])*
	 aluminum(wave[j])*atm(wave[j],1.3);
    a2 = bt*ip21(wave[j])*bk7(wave[j])*aluminum(wave[j])*aluminum(wave[j])
         *atm(wave[j],1.3);
    integral += 0.5*(y1+y2)*(wave[j]-wave[j-1]);
    A += 0.5*(a1+a2)*(wave[j]-wave[j-1]);
    a1 = a2;
    y1 = y2;
    j++;
  }

  return(integral/A);
}


double yflux(flux,wave)
double *flux,*wave;
{
  static double ywave[] = {5150.0,5175.0,5200.0,5225.0,5250.0,5275.0,5300.0,
     5325.0,5350.0,5375.0,5400.0,5425.0,5450.0,5475.0,5500.0,5525.0,5550.0,
     5575.0,5600.0,5625.0,5650.0,5675.0,5700.0,5725.0,5750.0,5775.0,5800.0,
     5825.0,5850.0};
  static double ytrans[] = {0.00,0.68,1.68,2.72,4.00,6.96,10.20,15.20,23.40,
	  33.00,41.00,45.40,47.92,50.12,52.20,51.32,45.44,34.80,24.80,15.48,
	  10.00,6.60,4.16,2.52,1.36,0.96,0.80,0.40,0.00};
  double Wave,y1,y2,yt,integral=0.0,A=0.0,a1,a2;
  int i;
  int k = 1;
  int j=0;


  y1 = 0.0;
  a1 = 0.0;
  j = 0;
  while(wave[j] < ywave[28]) {
    if(wave[j] <= ywave[0]) {
      j++;
      continue;
    }
    if(wave[j] >= ywave[28]) break;
    for(i=k;i<=28;i++) {
      if(wave[j] > ywave[i-1] && wave[j] <= ywave[i]) {
	k = i;
	break;
      }
    }
    yt =  ytrans[k-1] + (ytrans[k]-ytrans[k-1])*(wave[j]-ywave[k-1])/
	  (ywave[k]-ywave[k-1]);
    y2 = yt*flux[j]*ip21(wave[j])*bk7(wave[j])*aluminum(wave[j])*
	 aluminum(wave[j])*atm(wave[j],1.3);
    a2 = yt*ip21(wave[j])*bk7(wave[j])*aluminum(wave[j])*aluminum(wave[j])
         *atm(wave[j],1.3);
    integral += 0.5*(y1+y2)*(wave[j]-wave[j-1]);
    A += 0.5*(a1+a2)*(wave[j]-wave[j-1]);
    y1 = y2;
    a1 = a2;
    j++;
  }
  return(integral/A);
}

double ip21(wave)
double wave;
{
  static double wv[] = {3000.0,3100.0,3200.0,3300.0,3400.0,3500.0,3600.0,
   3700.0,3800.0,3900.0,4000.0,4100.0,4200.0,4300.0,4400.0,4500.0,4600.0,
   4700.0,4800.0,4900.0,5000.0,5100.0,5200.0,5300.0,5400.0,5500.0,5600.0,
   5700.0,5800.0,5900.0,6000.0};

  static double sens[] = 
   {0.08,0.22,0.45,0.68,0.81,0.90,0.95,0.97,0.99,1.00,1.00,
    0.99,0.98,0.95,0.91,0.87,0.83,0.77,0.71,0.65,0.58,0.52,0.46,0.40,0.34,
    0.29,0.24,0.20,0.16,0.13,0.10};
  int j,k;

  k = 0;
  for(j=0;j<31;j++) {
    if(wave <= wv[j]) {
      k = j;
      break;
    }
  }

  return(sens[k-1] + (sens[k]-sens[k-1])*(wave-wv[k-1])/(wv[k]-wv[k-1]));
}

double bk7(wave)
double wave;
{
  static double wv[] = 
  {3000.0,3100.0,3200.0,3300.0,3400.0,3500.0,3600.0,3700.0,
   3800.0,3900.0,4000.0,4200.0,4400.0,4600.0,4800.0,5000.0,5400.0,5800.0,
   6200.0,6600.0,7000.0};

  static double trans[] = 
 {0.26,0.59,0.81,0.93,0.97,0.986,0.992,0.995,0.996,0.998,
  0.998,0.998,0.999,0.999,0.999,0.999,0.999,0.999,0.999,0.999,0.999};
  int j,k;

  k = 0;
  for(j=0;j<21;j++) {
    if(wave <= wv[j]) {
      k = j;
      break;
    }
  }

  return(trans[k-1] + (trans[k]-trans[k-1])*(wave-wv[k-1])/(wv[k]-wv[k-1]));
}

double aluminum(wave)
double wave;
{
  static double wv[] = 
 {2952.0,3100.0,3263.0,3444.0,3646.0,3874.0,4133.0,4428.0,
  4768.0,5166.0,5635.0,6199.0,6525.0,6888.0};

  static double trans[] = {0.9248,0.9248,0.9247,0.9246,0.9245,0.9243,0.9241,
  0.9242,0.9238,0.9228,0.9200,0.9148,0.9116,0.9069};
  int j,k;

  k = 0;
  for(j=0;j<14;j++) {
    if(wave <= wv[j]) {
      k = j;
      break;
    }
  }

  return(trans[k-1] + (trans[k]-trans[k-1])*(wave-wv[k-1])/(wv[k]-wv[k-1]));
}

double atm(wave,a)
double wave;
double a;
{
  static double wa[] = 
                {2000.0,2200.0,2400.0,2600.0,2800.0,3000.0,3200.0,3400.0,
		 3600.0,3800.0,4000.0,4500.0,5000.0,5500.0,6000.0,
		 6500.0,7000.0};
  static double at[] = 
                {20.0,27.0,68.0,89.0,36.0,4.5,1.30,0.84,0.68,0.55,0.46,0.31,
		 0.23,0.195,0.170,0.126,0.092};
  int j,k;
  double ex;

  k = 0;
  for(j=0;j<17;j++) {
    if(wave <= wa[j]) {
      k = j;
      break;
    }
  }

  ex = at[k-1] + (at[k]-at[k-1])*(wave-wa[k-1])/(wa[k]-wa[k-1]);
  return(exp(-ex*a));
}