RiceMunk
commented
8 years ago Here's an ancient bit of C code that I'm planning on adapting for the KK relation:
/* From tielens Thu Dec 5 17:52:06 1996
Received: by dusty.arc.nasa.gov (8.6.8/1.35)
id RAA10340; Thu, 5 Dec 1996 17:52:06 -0800
From: tielens (A.G.G.M. Tielens)
Message-Id: <199612060152.RAA10340@dusty.arc.nasa.gov>
Subject: optcon.f revised
To: boogert (Adwin Boogert)
Date: Thu, 5 Dec 1996 17:52:06 -0800 (PST)
X-Mailer: ELM [version 2.4 PL22]
Content-Type: text
Content-Length: 3745
Status: RO
OPTCON.C: This program calculates optical constants
from transmission measurements using the Kramers-Kronig relationship.
v0.0: 27-12-1996 converted from Fortran to C (BOOGERT)
NOTE: The number of frequency points is defined by ARRDIM.
Adjust it according to length of input file.
Finally, note the sign of the KK integral ! */
#include<stdio.h>
#include<math.h>
#include"dcomplex.h"
#define ARRDIM 20000 /* dimension of arrays, must be < NPOINTS, i.e.
number of spectral datapoints */
#define PI 4.*atan(1.)
main()
{
dcomplex xc,yc,n1[ARRDIM],eix2,eix,r01,r12,t01,t12,t02,j ;
dcomplex m0, m1[ARRDIM], m2, Trans(), Refl();
double nu0,x,xx,alpha[ARRDIM],nu[ARRDIM],n01,d,delta ;
double t[ARRDIM],alpnew,tcal[ARRDIM], dif, zc, kkint(),dif1,dif2;
FILE *fpi, *fpo;
char filein[30];
char fileout[40];
int k, dk, nn,NPOINTS;
/* INPUT SECTION */
printf("Geef optische brekingsindex n01\n");
scanf("%lf", &n01 ); /* N01~1.22-1.35 op grote frequentie */
printf("Geef dikte ijs in micron\n");
scanf("%lf", &d ); /* lab ijsdikte D in cm, e.g. D=3.0D-5 */
printf("Geef input file naam\n");
scanf("%s", filein );
printf("Geef output file naam\n");
scanf("%s", fileout );
printf("Number of datapoints\n");
scanf("%i", &NPOINTS );
/*NUMBER OF DATAPOINTS MUST BE <= ARRAY DIMENSION*/
if (NPOINTS > ARRDIM) {
printf("\a\a\a **Array dimensions too small!!! \a\a\a\n");
printf("\a\a\a **Program aborted, adjust ARRDIM \a\a\a\n");
} else {
j=DComplex(0.,1.); /* the complex i */
m0=DComplex(1., 0.); /* complex refractive index of vacuum */
m2=DComplex(1.74, 0.); /* complex refractive index of CsI window*/
d=d*1.e-4; /*thickness ice sample in cm*/
/*printf("\t%f\t%f\n\n", n01,d);*/
delta=1.e-5 ; /*stop criterion final iteration, |1-Tcal/Tobs|*/
/*printf("\t%lf\n\n", delta);*/
/* READ DATA FROM FILE*/
fpi=fopen(filein,"r");
for (k=0; k < NPOINTS; k++) {
fscanf(fpi,"%le %le", &nu[k], &t[k]);
/*printf("%i got here\n\n\n", k);*/
/* printf("%e\t%e\n", nu[k],t[k]);*/
t[k]=pow(10,-t[k]);
m1[k]=DComplex(1.3,0.);
alpha[k]=0.;
}
nn=0;
dif=1.;
/* BEGINNING LOOP TO CONVERGE TO OBSERVED SPECTRUM */
while (dif > delta) {
nn++;
/* GIVEN N1(K), CALCULATE ALPHA (AND K1(K))
FROM MEASURED SPECTRA */
for (k=0; k<NPOINTS ; k++) {
xc=DCmul(DRCmul(4.*PI*d*nu[k], m1[k]), j);
eix=DCexp(xc);
r01=Refl(m0,m1[k]);
r12=Refl(m1[k],m2);
t01=Trans(m0,m1[k]);
t12=Trans(m1[k], m2);
t02=Trans(m0, m2);
yc=DCdiv(DCdiv(DCmul(t01,t12),t02),(DCadd(DComplex(1.,0.),DCmul(DCmul(r01,r12),eix))));
zc=DCmul(yc,DConjg(yc)).r; /* this is a double */
alpnew=(-log(t[k])+log(zc))/d;
alpha[k]=alpnew;
}
/* GIVEN ALPHA (OR K1(K)), M1(K) IS CALCULATED
USING THE KRAMERS KRONIG RELATIONSHIP.
CALCULATE SPECTRUM AND COMPARE WITH MEASURED ONE. */
dif=0.;
for (k=0; k<NPOINTS; k++) {
m1[k]=DComplex(kkint(k,nu, alpha, n01,NPOINTS),alpha[k]/(4.*PI*nu[k]));
if (dif <= delta) {
xc=DCmul(DRCmul(2.*PI*d*nu[k], m1[k]), j);
eix=DCexp(xc);
eix2=DCexp(DRCmul(2.,xc));
r01=Refl(m0,m1[k]);
r12=Refl(m1[k],m2);
t01=Trans(m0,m1[k]);
t12=Trans(m1[k], m2);
t02=Trans(m0, m2);
yc=DCdiv(DCdiv(DCmul(DCmul(t01,t12),eix),t02),(DCadd(DComplex(1.,0.),DCmul(DCmul(r01,r12),eix2))));
x=DCmul(yc,DConjg(yc)).r; /* this is a double */
tcal[k]=x;
if (x==0. || t[k]==0.)
dif2=fabs(x-t[k]);
else
dif2=fabs(1.-t[k]/x);
if (dif2 > dif) {
dif=dif2; /* the maximum deviating point ( dangerous!!)*/
dk=k;
}
/* printf("%le\t%le\n",dif2,dif); */
/* printf("%i\t",k);*/
}
}
printf("\a iteration %i difference %le at %le cm-1\n",nn,dif,nu[dk]);
} /* end of convergence while loop*/
printf("got here after %i iterations\n", nn);
fpo=fopen(fileout,"w");
for (k=0; k<NPOINTS; k++) {
fprintf(fpo, "%le\t%le\t%le\t%le\t%le\n", nu[k],m1[k].r,m1[k].i,t[k],tcal[k]);
}
}
}
/* FUNCTION "Refl": Complex reflection coefficient medium p-->q */
dcomplex Refl(Np, Nq)
dcomplex Np, Nq;
{
dcomplex r;
r=DCdiv(DCsub(Np,Nq),DCadd(Np,Nq));
return r;
}
/* FUNCTION "Trans": Complex transmission coefficient medium p-->q*/
dcomplex Trans(dcomplex Np, dcomplex Nq)
{
return DCdiv(DRCmul(2.,Np),DCadd(Np,Nq));
}
/* FUNCTION "KKINT": Kramers-Kroenig integration */
double kkint(k,nu, alpha, n01, NPOINTS)
double nu[ARRDIM], alpha[ARRDIM], n01 ;
int k, NPOINTS;
{
int i;
double dn, n, wnu[ARRDIM];
/* Note that sign in Kramers-Kronig equation depends
on sign of WNU (ie., deltanu=+ sign is +).
Use trapezium rule, but consider more sophisticated integration
routine in near future. */
for (i=1; i<(NPOINTS-1); i++) {
wnu[i]=(nu[i+1]-nu[i-1])/2.;
}
wnu[0]=(nu[1]-nu[0])/2.;
wnu[NPOINTS-1]=(nu[NPOINTS-1]-nu[NPOINTS-2])/2.;
dn=0.;
for (i=0; i<NPOINTS; i++) {
if (i != k)
dn=dn+alpha[i]/(pow(nu[i],2.)-pow(nu[k],2.))*wnu[i];
}
n=n01+dn/(2.*pow(PI,2));
/* printf("\t%.25lf\n\n", n);*/
return n ;
}
Currently the
CDEspectrumclass makes use of complex refractive index values which are commonly calculated using the Kramers-Kronig (KK) relation equation. It is not a terribly complex calculation, yet it's not included as part of Omnifit.An obvious improvement to Omnifit is to make it be included as part of the library, so
CDEspectrumcan be initialised without outside help.