023c8f6ae56090dcc2520df3cc44cdfc3e37a964
[physik/nlsop.git] / parse_trim_collision.c
1 /*
2  * parse srim 2003.26 collision data 
3  *
4  * author: hackbard@hackdaworld.dyndns.org
5  *
6  * this may just work with srim 2003.26 and 180keV c+ -> si,
7  * not quite sure though! :)
8  *
9  */
10
11 #define _GNU_SOURCE
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <unistd.h>
16 #include <sys/types.h>
17 #include <sys/stat.h>
18 #include <fcntl.h>
19
20 #define MAXION 5000
21 #define MAXZ 233
22
23 typedef struct s_z {
24         int cell_hits_from_ions;
25         int total_hits;
26         double energy;
27 } t_z;
28
29 int main(int argc,char **argv) {
30
31         int fd,ion,i,j;
32         int z,skipped;
33         int avg1,avg2;
34         double en;
35         double max_e;
36         int max_chfi,max_th;
37         t_z Z[MAXZ];
38         unsigned char hit[MAXZ];
39         char buf[256],*p,value[6],value2[11];
40
41         i=0;
42
43         if(argc!=2) return 0;
44         
45         printf("opening file %s ...\n",argv[1]);
46         if((fd=open(argv[1],O_RDONLY))<0) {
47                 printf("unable to open file %s\n",argv[1]);
48                 return 0;
49         }
50
51         ion=1;
52         skipped=0;
53         avg1=0;
54         avg2=0;
55         for(i=0;i<MAXZ;i++) {
56                 Z[i].cell_hits_from_ions=0;
57                 Z[i].total_hits=0;
58                 Z[i].energy=0;
59         }
60         memset(hit,0,MAXZ);
61
62         while(1) {
63
64                 /* read line */
65                 i=0;
66                 memset(buf,0,256);
67                 while(1) {
68                         j=read(fd,&buf[i],1);
69                         if(j<=0) {
70                                 close(fd);
71                                 /* norm */
72                                 max_e=Z[0].energy;
73                                 max_th=Z[0].total_hits;
74                                 max_chfi=Z[0].cell_hits_from_ions;
75                                 for(i=0;i<MAXZ;i++) {
76                                         if(Z[i].energy>max_e) max_e=Z[i].energy;
77                                         if(Z[i].total_hits>max_th) max_th=Z[i].total_hits;
78                                         if(Z[i].cell_hits_from_ions>max_chfi) max_chfi=Z[i].cell_hits_from_ions;
79                                         avg1+=Z[i].total_hits;
80                                         avg2+=Z[i].cell_hits_from_ions;
81                                 }
82                                 for(i=0;i<MAXZ;i++) 
83                                         //printf("%d %f %f %f\n",i*3,1.0*Z[i].total_hits/max_th,
84                                         //                      Z[i].energy/max_e,
85                                         //                      1.0*Z[i].cell_hits_from_ions/max_chfi);
86                                         printf("%d %f\n",i*3,1.0*Z[i].total_hits/max_th); // hits
87                                         //printf("%d %f\n",i*3,1.0*Z[i].energy/max_e); // energy
88                                 printf("skipped = %d\n",skipped);
89                                 printf("average hits per ion %d / %d\n",avg1/ion,avg2/ion);
90                                 return 1;
91                         }
92                         if(buf[i]=='\n') break;
93                         i++;
94                 }
95
96                 /* parse line */
97                 if((buf[0]=='³')&&(buf[1]!='=')&&(buf[1]!=' ')) {
98                         p=strtok(buf,"³");
99                         value[0]=p[0];
100                         value[1]=p[1];
101                         value[2]=p[2];
102                         value[3]=p[3];
103                         value[4]=p[4];
104                         value[5]='\0';
105                         p=strtok(NULL,"³");
106                         value2[0]=p[0];
107                         value2[1]=p[1];
108                         value2[2]='.';
109                         value2[3]=p[3];
110                         value2[4]=p[4];
111                         value2[5]=p[5];
112                         value2[6]=p[6];
113                         value2[7]=p[7];
114                         value2[8]=p[8];
115                         value2[9]='\0';
116                         en=atof(value2);
117                         p=strtok(NULL,"³");
118                         value2[0]=p[0];
119                         value2[1]=p[1];
120                         value2[2]=p[2];
121                         value2[3]=p[3];
122                         value2[4]=p[4];
123                         value2[5]='.';
124                         value2[6]=p[6];
125                         value2[7]=p[7];
126                         value2[8]=p[8];
127                         value2[9]=p[9];
128                         value2[10]='\0';
129                         z=(int)(atof(value2)/30.);
130                         if(z>232) skipped+=1;
131                         else {
132                                 //Z[z].energy+=en;
133                                 //Z[z].total_hits+=1;
134                                 hit[z]=1;
135                         }
136                         if(ion!=atoi(value)) {
137                                 /* new ion */
138                                 for(i=0;i<MAXZ;i++) Z[i].cell_hits_from_ions+=hit[i];
139                                 memset(hit,0,MAXZ);
140                                 ion=atoi(value);
141                         }
142                         // printf("%d %d %f %d\n",z*3,Z[z].total_hits,
143                         //      Z[z].energy,Z[z].cell_hits_from_ions);
144                         // return 0;
145                 }
146                 if((buf[0]=='Û')&&(buf[1]==' ')) {
147                         value2[0]=buf[25];
148                         value2[1]=buf[26];
149                         value2[2]=buf[27];
150                         value2[3]=buf[28];
151                         value2[4]='.';
152                         value2[5]=buf[30];
153                         value2[6]=buf[31];
154                         value2[7]=buf[32];
155                         value2[8]=buf[33];
156                         value2[9]='\0';
157                         z=(int)(atof(value2)/30.);
158                         value2[0]=buf[14];
159                         value2[1]=buf[15];
160                         value2[2]=buf[16];
161                         value2[3]=buf[17];
162                         value2[4]=buf[18];
163                         value2[5]='.';
164                         value2[6]=buf[20];
165                         value2[7]=buf[21];
166                         value2[8]=buf[22];
167                         value2[9]=buf[23];
168                         value2[10]='\0';
169                         en=atof(value2);
170                         if(z>232) skipped+=1;
171                         else {
172                                 Z[z].energy+=en;
173                                 Z[z].total_hits+=1;
174                                 hit[z]=1;
175                         }
176                 }
177         }
178         return 0;
179 }