##// END OF EJS Templates
test 4
rflores -
r1637:790a364b06d6
parent child
Show More
@@ -1,282 +1,282
1 #include <stdio.h>
1 #include <stdio.h>
2 #include <math.h>
2 #include <math.h>
3 #include <string.h>
3 #include <string.h>
4 #include <stdlib.h>
4 #include <stdlib.h>
5 #include <fcntl.h>
5 #include <fcntl.h>
6 #include <time.h>
6 #include <time.h>
7 #include "complex.h"
7 #include "complex.h"
8
8
9 #define NFREQ 512
9 #define NFREQ 512
10 #define LNES 12
10 #define LNES 12
11 #define LTES 16
11 #define LTES 16
12 #define NANG 19
12 #define NANG 19
13
13
14 #define MAX(x,y) ((x)>(y) ? (x) : (y))
14 #define MAX(x,y) ((x)>(y) ? (x) : (y))
15 #define MIN(x,y) ((x)<(y) ? (x) : (y))
15 #define MIN(x,y) ((x)<(y) ? (x) : (y))
16
16
17 fcomplex ****exlib;
17 fcomplex ****exlib;
18 int first=1;
18 int first=1;
19
19
20 void swab_test(void);
20 void swab_test(void);
21 //float swab(signed char *a);
21 //float swab(signed char *a);
22 void read_exlib(char *s,int nfreq, int lnes, int lTes, int nang,
22 void read_exlib(char *s,int nfreq, int lnes, int lTes, int nang,
23 fcomplex ****exlib);
23 fcomplex ****exlib);
24
24
25 void initialize(void);
25 void initialize(void);
26 //void collision_(float *np, float *tp, float *fp, float *ap, float *yr, float *yi)// step 1
26 //void collision_(float *np, float *tp, float *fp, float *ap, float *yr, float *yi)// step 1
27 void collision_(float *n, float *t, float *f, float *a, fcomplex *y);
27 void collision_(float *n, float *t, float *f, float *a, fcomplex *y);
28
28
29 /*
29 /*
30 void main(void){
30 void main(void){
31 float n=4.9e11,t=1100.0,f=300.0,a=1.0;
31 float n=4.9e11,t=1100.0,f=300.0,a=1.0;
32 int i;
32 int i;
33 fcomplex j;
33 fcomplex j;
34 for(i=0;i<10;i++){
34 for(i=0;i<10;i++){
35 a+=0.05;
35 a+=0.05;
36 collision_(&n,&t,&f,&a,&j);
36 collision_(&n,&t,&f,&a,&j);
37 printf("%f %f %f \n",a,j.r,j.i);
37 printf("%f %f %f \n",a,j.r,j.i);
38 }
38 }
39 }
39 }
40 */
40 */
41
41
42 //void collision_(float *np, float *tp, float *fp, float *ap, float *yr, float *yi){ // step 2
42 //void collision_(float *np, float *tp, float *fp, float *ap, float *yr, float *yi){ // step 2
43 void collision_(float *np, float *tp, float *fp, float *ap, fcomplex *y){
43 void collision_(float *np, float *tp, float *fp, float *ap, fcomplex *y){
44 float ne[]={5.0e10, 1.0e11, 2.0e11, 3.5e11, 5.0e11, 7.5e11,
44 float ne[]={5.0e10, 1.0e11, 2.0e11, 3.5e11, 5.0e11, 7.5e11,
45 1.0e12, 1.5e12, 2.0e12, 2.5e12, 3.0e12, 3.5e12};
45 1.0e12, 1.5e12, 2.0e12, 2.5e12, 3.0e12, 3.5e12};
46 float te[]={600.0, 700.0, 800.0, 900.0, 1000.0, 1150.0, 1300.0, 1500.0,
46 float te[]={600.0, 700.0, 800.0, 900.0, 1000.0, 1150.0, 1300.0, 1500.0,
47 1750.0, 2000.0, 2250.0, 2500.0, 2750.0, 3000.0, 3500.0, 4000.0};
47 1750.0, 2000.0, 2250.0, 2500.0, 2750.0, 3000.0, 3500.0, 4000.0};
48 float alpha[]={.125, .1875, .25, .3125, .375, .5, .625, .75, 1.0, 1.25,
48 float alpha[]={.125, .1875, .25, .3125, .375, .5, .625, .75, 1.0, 1.25,
49 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0};
49 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0};
50 float freq[2],fmax=6237.79;
50 float freq[2],fmax=6237.79;
51 float n,t,f,a;
51 float n,t,f,a;
52 float v,w1,w2,w3,w4,weight;
52 float v,w1,w2,w3,w4,weight;
53 fcomplex g;
53 fcomplex g;
54 int i,j,k,l,m,i1,j1,k1,l1;
54 int i,j,k,l,m,i1,j1,k1,l1;
55
55
56 if(first==1){
56 if(first==1){
57 //printf("initi");
57 //printf("initi");
58 initialize();
58 initialize();
59 first=0;
59 first=0;
60 }
60 }
61 //printf(" inputs in collision_ \n y.r %f y.i %f \n *************\n",(*y).r,(*y).i);
61 //printf(" inputs in collision_ \n y.r %f y.i %f \n *************\n",(*y).r,(*y).i);
62
62
63 //printf("After:print\n");
63 //printf("After:print\n");
64 a=*ap; t=*tp; f=*fp; n=*np;
64 a=*ap; t=*tp; f=*fp; n=*np;
65
65
66 a=MAX(MIN(a,6.0),0.125);
66 a=MAX(MIN(a,6.0),0.125);
67 t=MAX(MIN(t,4000.0),600.0);
67 t=MAX(MIN(t,4000.0),600.0);
68 f=MIN(f,fmax-1.0);
68 f=MIN(f,fmax-1.0);
69 n=MAX(MIN(n,3.5e12),5.0e10);
69 n=MAX(MIN(n,3.5e12),5.0e10);
70
70
71 i1=((float)(NFREQ-1))*f/fmax;
71 i1=((float)(NFREQ-1))*f/fmax;
72 i1=MIN(i1,NFREQ-1);
72 i1=MIN(i1,NFREQ-1);
73 freq[0]=i1*fmax/(float)(NFREQ-1);
73 freq[0]=i1*fmax/(float)(NFREQ-1);
74 freq[1]=(i1+1)*fmax/(float)(NFREQ-1);
74 freq[1]=(i1+1)*fmax/(float)(NFREQ-1);
75
75
76 j1=0;
76 j1=0;
77 while(n>ne[j1+1]) j1++;
77 while(n>ne[j1+1]) j1++;
78 j1=MIN(j1,LNES-1);
78 j1=MIN(j1,LNES-1);
79
79
80 k1=0;
80 k1=0;
81 while(t>te[k1+1]) k1++;
81 while(t>te[k1+1]) k1++;
82 k1=MIN(k1,LTES);
82 k1=MIN(k1,LTES);
83
83
84 l1=0;
84 l1=0;
85 while(a>alpha[l1+1]) l1++;
85 while(a>alpha[l1+1]) l1++;
86 l1=MIN(l1,NANG);
86 l1=MIN(l1,NANG);
87
87
88 // printf("i,j,k,l: %i %i %i %i \n",i1,j1,k1,l1);
88 // printf("i,j,k,l: %i %i %i %i \n",i1,j1,k1,l1);
89
89
90 v=(fmax/(float)(NFREQ-1))*
90 v=(fmax/(float)(NFREQ-1))*
91 (ne[j1+1]-ne[j1])*(te[k1+1]-te[k1])*(alpha[l1+1]-alpha[l1]);
91 (ne[j1+1]-ne[j1])*(te[k1+1]-te[k1])*(alpha[l1+1]-alpha[l1]);
92 g=Complex(0.0,0.0);
92 g=Complex(0.0,0.0);
93
93
94 for(i=0;i<2;i++){
94 for(i=0;i<2;i++){
95 if(i==0)
95 if(i==0)
96 w1=freq[1]-f;
96 w1=freq[1]-f;
97 else
97 else
98 w1=f-freq[0];
98 w1=f-freq[0];
99 for(j=0;j<2;j++){
99 for(j=0;j<2;j++){
100 if(j==0)
100 if(j==0)
101 w2=ne[j1+1]-n;
101 w2=ne[j1+1]-n;
102 else
102 else
103 w2=n-ne[j1];
103 w2=n-ne[j1];
104 for(k=0;k<2;k++){
104 for(k=0;k<2;k++){
105 if(k==0)
105 if(k==0)
106 w3=te[k1+1]-t;
106 w3=te[k1+1]-t;
107 else
107 else
108 w3=t-te[k1];
108 w3=t-te[k1];
109 for(l=0;l<2;l++){
109 for(l=0;l<2;l++){
110 if(l==0)
110 if(l==0)
111 w4=alpha[l1+1]-a;
111 w4=alpha[l1+1]-a;
112 else
112 else
113 w4=a-alpha[l1];
113 w4=a-alpha[l1];
114 weight=w1*w2*w3*w4;
114 weight=w1*w2*w3*w4;
115 g.r+=exlib[i1+i][j1+j][k1+k][l1+l].r*weight;
115 g.r+=exlib[i1+i][j1+j][k1+k][l1+l].r*weight;
116 g.i+=exlib[i1+i][j1+j][k1+k][l1+l].i*weight;
116 g.i+=exlib[i1+i][j1+j][k1+k][l1+l].i*weight;
117 // printf(" in collision i1+i %d j1+j %d k1+k %d l1+l %d exlib.r %f exlib.i %f \n",i1+i,j1+j,k1+k,l1+l,exlib[i1+i][j1+j][k1+k][l1+l].r,exlib[i1+i][j1+j][k1+k][l1+l].i);
117 // printf(" in collision i1+i %d j1+j %d k1+k %d l1+l %d exlib.r %f exlib.i %f \n",i1+i,j1+j,k1+k,l1+l,exlib[i1+i][j1+j][k1+k][l1+l].r,exlib[i1+i][j1+j][k1+k][l1+l].i);
118 // printf(" in collision g.r %f g.i %f \n",g.r,g.i);
118 // printf(" in collision g.r %f g.i %f \n",g.r,g.i);
119 }
119 }
120 }
120 }
121 }
121 }
122 }
122 }
123 g.r/=v; g.i/=v;
123 g.r/=v; g.i/=v;
124
124
125 *y=g;
125 *y=g;
126 //*yr=g.r
126 //*yr=g.r
127 //*yi=g.i
127 //*yi=g.i
128 //printf(" outputs in collision_ \n g.r %f g.i %f \n *************\n",g.r,g.i);
128 //printf(" outputs in collision_ \n g.r %f g.i %f \n *************\n",g.r,g.i);
129 //exit(-2);
129 //exit(-2);
130 //getchar();
130 //getchar();
131 return;
131 return;
132 }
132 }
133
133
134 void initialize(void){
134 void initialize(void){
135 // char s[]="/usr/local/lib/faraday/jlib26feb2002";
135 // char s[]="/usr/local/lib/faraday/jlib26feb2002";
136 // char s[]="/usr/local/lib/faraday/jlib26feb2002";
136 // char s[]="/usr/local/lib/faraday/jlib26feb2002";
137 // FILE *s=fopen("/usr/local/lib/faraday/jlib26feb2002","r");
137 // FILE *s=fopen("/usr/local/lib/faraday/jlib26feb2002","r");
138 int nfreq=NFREQ, lnes=LNES, lTes=LTES, nang=NANG;
138 int nfreq=NFREQ, lnes=LNES, lTes=LTES, nang=NANG;
139 int i,j,k;
139 int i,j,k;
140 int crr,crr2;
140 int crr,crr2;
141 int the_number=1025;
141 int the_number=1025;
142
142
143 char the_path[1025];
143 char the_path[1025];
144
144
145 get_path_reader_(the_path,&the_number);
145 get_path_reader_(the_path,&the_number);
146 //printf("C the_number: %i",the_number);
146 //printf("C the_number: %i",the_number);
147 //printf("After\n");
147 //printf("After\n");
148 //printf("The Path: %s", the_path);
148 //printf("The Path: %s", the_path);
149 //printf("END\n");
149 //printf("END\n");
150 //printf("%s%n",the_path,&crr);
150 //printf("%s%n",the_path,&crr);
151 //printf("=%d",crr);
151 //printf("=%d",crr);
152 //printf("END2\n");
152 //printf("END2\n");
153 char the_path_true[the_number];
153 char the_path_true[the_number];
154 memcpy(the_path_true,the_path,the_number);
154 memcpy(the_path_true,the_path,the_number);
155 the_path_true[the_number] = '\0';
155 the_path_true[the_number] = '\0';
156 //printf("%s%n",the_path_true,&crr2);
156 printf("%s%n",the_path_true,&crr2);
157 //printf("=%d",crr2);
157 //printf("=%d",crr2);
158 //printf("END3\n");
158 //printf("END3\n");
159 //printf("PATH: %s",s);
159 //printf("PATH: %s",s);
160 //fprintf(s);
160 //fprintf(s);
161 exlib=malloc(nfreq*8);
161 exlib=malloc(nfreq*8);
162 for(i=0;i<nfreq;i++){
162 for(i=0;i<nfreq;i++){
163 exlib[i]=malloc(lnes*8);
163 exlib[i]=malloc(lnes*8);
164 for(j=0;j<lnes;j++){
164 for(j=0;j<lnes;j++){
165 exlib[i][j]=malloc(lTes*8);
165 exlib[i][j]=malloc(lTes*8);
166 for(k=0;k<lTes;k++){
166 for(k=0;k<lTes;k++){
167 exlib[i][j][k]=(fcomplex *)malloc(nang*sizeof(fcomplex));
167 exlib[i][j][k]=(fcomplex *)malloc(nang*sizeof(fcomplex));
168 }
168 }
169 }
169 }
170 }
170 }
171
171
172 read_exlib(the_path_true,nfreq,lnes,lTes,nang,exlib);
172 read_exlib(the_path_true,nfreq,lnes,lTes,nang,exlib);
173 }
173 }
174
174
175 float swab(signed char *a)
175 float swab(signed char *a)
176 // table is in big endian format
176 // table is in big endian format
177 {
177 {
178 float mant;
178 float mant;
179 unsigned char *amant;
179 unsigned char *amant;
180 printf(" in swab 0 *** a[0] %c a[1] %c a[2] %c a[3] %c\n",a[0],a[1],a[2],a[3]);
180 printf(" in swab 0 *** a[0] %c a[1] %c a[2] %c a[3] %c\n",a[0],a[1],a[2],a[3]);
181 amant=(unsigned char *) &mant;
181 amant=(unsigned char *) &mant;
182 amant[3]=a[0]; amant[2]=a[1]; amant[1]=a[2]; amant[0]=a[3];
182 amant[3]=a[0]; amant[2]=a[1]; amant[1]=a[2]; amant[0]=a[3];
183 //printf(" in swab 1 ***\n");
183 //printf(" in swab 1 ***\n");
184 return(mant);
184 return(mant);
185 }
185 }
186
186
187 void swab_test(void)
187 void swab_test(void)
188 // table is in big endian format
188 // table is in big endian format
189 {
189 {
190 printf("*********in swab stest ********\n");
190 printf("*********in swab stest ********\n");
191 getchar();
191 getchar();
192 }
192 }
193
193
194
194
195 void read_exlib(char *the_path_true,int nfreq, int lnes, int lTes, int nang,
195 void read_exlib(char *the_path_true,int nfreq, int lnes, int lTes, int nang,
196 fcomplex ****exlib)
196 fcomplex ****exlib)
197 {
197 {
198 FILE *fp;
198 FILE *fp;
199 int i,j,k,l;
199 int i,j,k,l;
200 float a;
200 float a;
201 float t1, t2;
201 float t1, t2;
202 signed char *t1a,*t2b;
202 signed char *t1a,*t2b;
203 float mant;
203 float mant;
204 unsigned char *amant;
204 unsigned char *amant;
205 printf("Start reading exlib file \n");
205 printf("Start reading exlib file \n");
206 if( (fp=fopen(the_path_true,"r")) == NULL){
206 if( (fp=fopen(the_path_true,"r")) == NULL){
207 printf("Je Library file opening error\n");
207 printf("Je Library file opening error\n");
208 //printf("%s\n", strerror(errno));
208 //printf("%s\n", strerror(errno));
209 exit(-2);
209 exit(-2);
210 }
210 }
211 for(l=0;l<nang;l++)
211 for(l=0;l<nang;l++)
212 for(k=0;k<lTes;k++)
212 for(k=0;k<lTes;k++)
213 for(j=0;j<lnes;j++)
213 for(j=0;j<lnes;j++)
214 for(i=0;i<nfreq;i++){
214 for(i=0;i<nfreq;i++){
215 //printf(" in read_exlib 0***\n");
215 //printf(" in read_exlib 0***\n");
216 fread(&t1,sizeof(float),1,fp);
216 fread(&t1,sizeof(float),1,fp);
217 fread(&t2,sizeof(float),1,fp);
217 fread(&t2,sizeof(float),1,fp);
218 //printf(" in read_exlib 1***\n");
218 //printf(" in read_exlib 1***\n");
219 //printf(" i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,t1,t2);
219 //printf(" i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,t1,t2);
220 t1a=(signed char *)&t1;
220 t1a=(signed char *)&t1;
221 t2b=(signed char *)&t2;
221 t2b=(signed char *)&t2;
222 //printf(" in read_exlib 1.5 t1a %d t2b %d ***\n",t1a,t2b);
222 //printf(" in read_exlib 1.5 t1a %d t2b %d ***\n",t1a,t2b);
223 ////t1=swab(t1a);
223 ////t1=swab(t1a);
224 //printf(" t1 -- in readexlib 1.8 0 *** t1a[0] %c t1a[1] %c t1a[2] %c t1a[3] %c\n",t1a[0],t1a[1],t1a[2],t1a[3]);
224 //printf(" t1 -- in readexlib 1.8 0 *** t1a[0] %c t1a[1] %c t1a[2] %c t1a[3] %c\n",t1a[0],t1a[1],t1a[2],t1a[3]);
225 amant=(unsigned char *) &mant;
225 amant=(unsigned char *) &mant;
226 amant[3]=t1a[0]; amant[2]=t1a[1]; amant[1]=t1a[2]; amant[0]=t1a[3];
226 amant[3]=t1a[0]; amant[2]=t1a[1]; amant[1]=t1a[2]; amant[0]=t1a[3];
227 t1=mant;
227 t1=mant;
228
228
229 ////t2=swab(t2b);
229 ////t2=swab(t2b);
230 //printf(" t2 -- in readexlib 1.8 0 *** t2b[0] %c t2b[1] %c t2b[2] %c t2b[3] %c\n",t2b[0],t2b[1],t2b[2],t2b[3]);
230 //printf(" t2 -- in readexlib 1.8 0 *** t2b[0] %c t2b[1] %c t2b[2] %c t2b[3] %c\n",t2b[0],t2b[1],t2b[2],t2b[3]);
231 amant=(unsigned char *) &mant;
231 amant=(unsigned char *) &mant;
232 amant[3]=t2b[0]; amant[2]=t2b[1]; amant[1]=t2b[2]; amant[0]=t2b[3];
232 amant[3]=t2b[0]; amant[2]=t2b[1]; amant[1]=t2b[2]; amant[0]=t2b[3];
233 t2=mant;
233 t2=mant;
234
234
235 ////t1=swab((signed char *)&t1);
235 ////t1=swab((signed char *)&t1);
236 ////t2=swab((signed char *)&t2);
236 ////t2=swab((signed char *)&t2);
237 //printf(" in read_exlib 2 i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,t1,t2);
237 //printf(" in read_exlib 2 i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,t1,t2);
238 exlib[i][j][k][l]=Complex(t1,t2);
238 exlib[i][j][k][l]=Complex(t1,t2);
239 //printf(" in read_exlib 3 i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,exlib[i][j][k][l].r,exlib[i][j][k][l].i);
239 //printf(" in read_exlib 3 i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,exlib[i][j][k][l].r,exlib[i][j][k][l].i);
240 //getchar();
240 //getchar();
241 }
241 }
242 // getchar();
242 // getchar();
243
243
244 fclose(fp);
244 fclose(fp);
245 printf("Done reading exlib file \n");
245 printf("Done reading exlib file \n");
246 }
246 }
247 /*
247 /*
248 void read_exlib(char *s,int nfreq, int lnes, int lTes, int nang,
248 void read_exlib(char *s,int nfreq, int lnes, int lTes, int nang,
249 fcomplex ****exlib)
249 fcomplex ****exlib)
250 {
250 {
251 FILE *fp;
251 FILE *fp;
252 int i,j,k,l;
252 int i,j,k,l;
253 float a;
253 float a;
254 float t1, t2;
254 float t1, t2;
255
255
256 printf("Start reading exlib fileee \n");
256 printf("Start reading exlib fileee \n");
257
257
258 if( (fp=fopen(s,"r")) == NULL){
258 if( (fp=fopen(s,"r")) == NULL){
259 printf("Je Library file opening error\n");
259 printf("Je Library file opening error\n");
260 exit(-2);
260 exit(-2);
261 }
261 }
262
262
263 for(l=0;l<nang;l++)
263 for(l=0;l<nang;l++)
264 for(k=0;k<lTes;k++)
264 for(k=0;k<lTes;k++)
265 for(j=0;j<lnes;j++)
265 for(j=0;j<lnes;j++)
266 for(i=0;i<nfreq;i++){
266 for(i=0;i<nfreq;i++){
267 printf(" in read_exlib 0***\n");
267 printf(" in read_exlib 0***\n");
268 fread(&t1,sizeof(float),1,fp);
268 fread(&t1,sizeof(float),1,fp);
269 fread(&t2,sizeof(float),1,fp);
269 fread(&t2,sizeof(float),1,fp);
270 printf(" in read_exlib 1***\n");
270 printf(" in read_exlib 1***\n");
271 printf(" i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,t1,t2);
271 printf(" i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,t1,t2);
272 t1=swab((char *)&t1); t2=swab((char *)&t2);
272 t1=swab((char *)&t1); t2=swab((char *)&t2);
273 printf(" in read_exlib 2 i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,t1,t2);
273 printf(" in read_exlib 2 i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,t1,t2);
274 exlib[i][j][k][l]=Complex(t1,t2);
274 exlib[i][j][k][l]=Complex(t1,t2);
275 printf(" in read_exlib 3 i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,exlib[i][j][k][l].r,exlib[i][j][k][l].i);
275 printf(" in read_exlib 3 i %d j %d k %d l %d t1 %f t2 %f \n",i,j,k,l,exlib[i][j][k][l].r,exlib[i][j][k][l].i);
276 getchar();
276 getchar();
277 }
277 }
278
278
279 fclose(fp);
279 fclose(fp);
280 printf("Done reading exlib file \n");
280 printf("Done reading exlib file \n");
281 }
281 }
282 */
282 */
General Comments 0
You need to be logged in to leave comments. Login now