1 #include "mrilib.h"
2
3 #define sinc(x) (((x)==0.0) ? 1.0 : sin(x)/(x))
4
main(int argc,char * argv[])5 int main( int argc , char *argv[] )
6 {
7 int iarg , nP , nN , pp,nn ;
8 double rr,qq, lam , val , beta=0.0 ;
9
10 /* help? */
11
12 if( argc < 4 || strcmp(argv[1],"-help") == 0 ){
13 printf("Usage: sink N P lam [beta]\n"
14 "Outputs to stdout the .1D file with 2*P+1 columns and 2*N+1 rows:\n"
15 " sinc[ pi * (p-P*lam*n/N) ]\n"
16 "OR, if beta != 0:\n"
17 " sinc[ pi * (p-P*lam*( beta+(1-beta)*n/N)) ]\n"
18 " +sinc[ pi * (p-P*lam*(-beta+(1-beta)*n/N)) ]\n"
19 "for p=-P..P (column index) and n=-N..N (row index).\n"
20 ) ;
21 exit(0) ;
22 }
23
24 nN = (int)strtod(argv[1],NULL); if( nN < 1 ){fprintf(stderr,"illegal N\n");exit(1);}
25 nP = (int)strtod(argv[2],NULL); if( nP < 1 ){fprintf(stderr,"illegal P\n");exit(1);}
26 lam = strtod(argv[3],NULL);
27 if( argc > 4 ) beta = strtod(argv[4],NULL) ;
28
29 rr = (nP*lam*(1.0-beta))/nN ;
30 qq = (nP*lam*beta) ;
31
32 for( nn=-nN ; nn <= nN ; nn++ ){
33 for( pp=-nP ; pp <= nP ; pp++ ){
34 val = sinc( PI * (pp-rr*nn-qq) ) ;
35 #if 1
36 if( beta != 0.0 ) val += 0.737*sinc( PI * (pp-rr*nn+qq) ) ;
37 #endif
38 printf(" %.12g",val) ;
39 }
40 printf("\n") ;
41 }
42
43 exit(0) ;
44 }
45