Zcosft2.cin
// zcosft2.cin is the include file for compilation of the examples with the discrete cosine transforms, namely, DCTII and DCTIII.
// The first argument is array of length for some natural number $q$. If the numeration of elements of the array begins with zero, and the name of array is A, then, "A-1" should appear as first argument. (As usulally in C++, the zero'th argument is name of function.) The z_type should be defined as "double" in order to deal with real numbers or complex(double) in order to deal with complex input functions. FOr the applications in paraxial optics, the second case may be more suitable.
//The Second argument is just number of points of the array, it should be mentioned above.
// The Third (and the last) argument should be "-1" to perform DCTII and "1" to perform DCTIII.
// The results of the calculation are stored in the same array, specified as first argument.
void zcosft2(z_type y[], int n, int isign)
{
void zrealft(z_type data[], unsigned long n, int isign);
int i;
z_type sum,sum1,y1,y2,ytemp;
double theta,wi=0.0,wi1,wpi,wpr,wr=1.0,wr1,wtemp;
theta=0.5*M_PI/n;
wr1=cos(theta);
wi1=sin(theta);
wpr = -2.0*wi1*wi1;
wpi=sin(2.0*theta);
if(isign == 1){
for (i=1;i<=n/2;i++){
y1=0.5*(y[i]+y[n-i+1]);
y2=wi1*(y[i]-y[n-i+1]);
y[i]=y1+y2;
y[n-i+1]=y1-y2;
wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1;
wi1=wi1*wpr+wtemp*wpi+wi1;
}
zrealft(y,n,1);
for (i=3;i<=n;i+=2) {
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
y1=y[i]*wr-y[i+1]*wi;
y2=y[i+1]*wr+y[i]*wi;
y[i]=y1;
y[i+1]=y2;
}
sum=0.5*y[2];
for (i=n;i>=2;i-=2) {
sum1=sum;
sum += y[i];
y[i]=sum1;
}
} else if (isign == -1) {
ytemp=y[n];
for (i=n;i>=4;i-=2) y[i]=y[i-2]-y[i];
y[2]=2.0*ytemp;
for (i=3;i<=n;i+=2) {
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
y1=y[i]*wr+y[i+1]*wi;
y2=y[i+1]*wr-y[i]*wi;
y[i]=y1;
y[i+1]=y2;
}
zrealft(y,n,-1);
for(i=1;i<=n/2;i++) {
y1=y[i]+y[n-i+1];
y2=(0.5/wi1)*(y[i]-y[n-i+1]);
y[i]=0.5*(y1+y2);
y[n-i+1]=0.5*(y1-y2);
wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1;
wi1=wi1*wpr+wtemp*wpi+wi1;
}
}
}