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