git.gag.com Git - debian/atlc/blob - src/non_gui/do_fd_calculation.c

   1 /* atlc - arbitrary transmission line calculator, for the analysis of
   2 extern int number_of_workers;
   3 transmission lines are directional couplers.
   4
   5 Copyright (C) 2002. Dr. David Kirkby, PhD (G8WRB).
   6
   7 This program is free software; you can redistribute it and/or
   8 modify it under the terms of the GNU General Public License
   9 as published by the Free Software Foundation; either package_version 2
  10 of the License, or (at your option) any later package_version.
  11
  12 This program is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with this program; if not, write to the Free Software
  19 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
  20 USA.
  21
  22 Dr. David Kirkby, e-mail drkirkby at ntlworld.com
  23
  24 */
  25 #include "config.h"
  26
  27
  28 #ifdef HAVE_STDLIB_H
  29 #include <stdlib.h>
  30 #endif
  31
  32 #include "definitions.h"
  33
  34 extern int append_flag;
  35 extern int dielectrics_to_consider_just_now, coupler;
  36 extern int num_pes;
  37 extern double **Vij;
  38 extern int height;
  39 extern int number_of_workers;
  40
  41 void do_fd_calculation(struct transmission_line_properties *data, size_t size, FILE *where_to_print_fp, char *inputfile_filename)
  42 {
  43   double capacitance_old, capacitance;
  44   double velocity_of_light_in_vacuum;
  45   int count=0;
  46   /*
  47   if (data->dielectrics_in_bitmap > 1) {
  48     fprintf(stderr,"\nSorry, but on the 15th October 2003 I was advised there is an\n");
  49     fprintf(stderr,"error in atlc when computing systems with multiple dielectrics.\n");
  50     fprintf(stderr,"So until this problem is fixed, the facility has been disabled.\n\n");
  51     fprintf(stderr,"I hope to release a new version shortly without this problem.\n");
  52     exit(1);
  53   }
  54   */
  55   /* The line of best fit of non_metalic_pixels vs iterations required
  56   is y=0.0011 * non_metallic_pixels  + 283.
  57
  58   I'll ensure finite_difference is called about 10x by using
  59   0.00011 * the number of non-metallic elements  +28 as the number
  60   of times finite_difference is called each time */
  61
  62
  63
  64   /* The following 10 lines are for a single dielectric 2 conductor line */
  65   if (data->couplerQ==FALSE)
  66   {
  67     if(data->verbose_level >= 2)
  68       printf("Solving assuming a vacuum dielectric\n");
  69     capacitance=VERY_LARGE; /* Can be anything large */
  70         dielectrics_to_consider_just_now=1;
  71         data->dielectrics_to_consider_just_now=1;
  72
  73     do /* Start a finite calculation */
  74     {
  75       capacitance_old=capacitance;
  76
  77 #ifdef ENABLE_POSIX_THREADS
  78       if (number_of_workers == 0)
  79         capacitance=finite_difference_single_threaded();
  80       else
  81         capacitance=finite_difference_multi_threaded();
  82 #else
  83       capacitance=finite_difference_single_threaded();
  84 #endif
  85
  86       data->C_vacuum=capacitance;
  87       data->C=capacitance;
  88       data->L_vacuum=MU_0*EPSILON_0/capacitance; /* Same as L in *ALL* cases */
  89       data->Zo_vacuum=sqrt(data->L_vacuum/data->C_vacuum);  /* Standard formaul for Zo */
  90       data->C=capacitance;
  91       if (data->dielectrics_in_bitmap == 1) /* Just get C by simple scaling of Er */
  92       {
  93         data->C=capacitance*data->found_this_dielectric;  /* Scaled by the single dielectric constant */
  94                 data->Er=data->found_this_dielectric;
  95       }
  96       else
  97                 data->Er=1.0;
  98       data->Zo=sqrt(data->L_vacuum/data->C);  /* Standard formula for Zo */
  99       data->Zodd=sqrt(data->L_vacuum/data->C);  /* Standard formula for Zo */
 100       velocity_of_light_in_vacuum=1.0/(sqrt(MU_0 * EPSILON_0)); /* around 3x10^8 m/s */
 101       data->velocity=1.0/pow(data->L_vacuum*data->C,0.5);
 102       data->velocity_factor=data->velocity/velocity_of_light_in_vacuum;
 103       data->relative_permittivity=sqrt(data->velocity_factor); /* ??? XXXXXX */
 104       if(data->verbose_level > 0 ) /* Only needed if intermediate results wanted.  */
 105         print_data_for_two_conductor_lines(*data, where_to_print_fp, inputfile_filename);
 106       count++;
 107     } while (fabs((capacitance_old-capacitance)/capacitance_old) > data->cutoff); /* end of FD loop */
 108     if(data->verbose_level >=4)
 109       printf("Total of %d iterations ( %d calls to finite_difference() )\n",ITERATIONS*count,count);
 110
 111     if((data->write_binary_field_imagesQ == TRUE || data->write_bitmap_field_imagesQ == TRUE) && data->dielectrics_in_bitmap==1 )
 112       write_fields_for_two_conductor_lines(inputfile_filename, *data, size);
 113     if(data->verbose_level == 0 && data->dielectrics_in_bitmap==1 )
 114       print_data_for_two_conductor_lines(*data, where_to_print_fp, inputfile_filename);
 115
 116     if ( data->dielectrics_in_bitmap >1)
 117     {
 118       /* We know the capacitance and inductance for the air spaced line
 119       as we calculated it above. Howerver, whilst the inductance
 120       is independant of the dielectric, the capacitance is not, so this
 121       has to be recalculated, taking care not to alter the inductance
 122       at all */
 123       if(data->verbose_level >= 2)
 124         printf("Now taking into account the permittivities of the different dielectrics for 2 conductors.\n");
 125
 126       dielectrics_to_consider_just_now=3; /* Any number > 1 */
 127       data->dielectrics_to_consider_just_now=2; /* Any number > 1 */
 128
 129       capacitance=VERY_LARGE; /* Can be anything large */
 130
 131       do /* Start a finite calculation */
 132       {
 133         capacitance_old=capacitance;
 134 #ifdef ENABLE_POSIX_THREADS
 135       if (number_of_workers == 0)
 136          capacitance=finite_difference_single_threaded();
 137       else
 138         capacitance=finite_difference_multi_threaded();
 139 #else
 140       capacitance=finite_difference_single_threaded();
 141 #endif
 142         data->C=capacitance;
 143         data->C_non_vacuum=capacitance;
 144         data->Zo=sqrt(data->L_vacuum/data->C_non_vacuum);  /* Standard formula for Zo */
 145         data->velocity=1.0/pow(data->L_vacuum*data->C_non_vacuum,0.5);
 146         data->velocity_factor=data->velocity/velocity_of_light_in_vacuum;
 147         data->relative_permittivity=sqrt(data->velocity_factor); /* ??? XXXXXX */
 148                 data->Er=data->C/data->C_vacuum;
 149         if(data->verbose_level > 0 ) /* Only needed if intermediate results wanted. */
 150           print_data_for_two_conductor_lines(*data, where_to_print_fp, inputfile_filename);
 151       } while (fabs((capacitance_old-capacitance)/capacitance_old) > data->cutoff); /* end of FD loop */
 152
 153       /* We must print the results now, but only bother if the verbose level was
 154       not not incremented on the command line, otherwide there will be two duplicate
 155       lines */
 156
 157       if (data->verbose_level == 0)
 158         print_data_for_two_conductor_lines(*data, where_to_print_fp, inputfile_filename);
 159       if(data->write_binary_field_imagesQ == TRUE || data->write_bitmap_field_imagesQ == TRUE)
 160         write_fields_for_two_conductor_lines(inputfile_filename, *data, size);
 161     }
 162   }
 163   else if (data->couplerQ==TRUE)
 164   {
 165     /* The properties of a couplers will be computed in 2 or 4 stages
 166     1) Compute the odd-mode impedance, assuming a vacuum dielectric, or
 167     if there is just one dielectric, that one.
 168
 169     2) Compute the odd-mode impedance, taking into account the effect of
 170     multiple dielectrics, IF NECESSARY
 171
 172     at this point, the negative voltages will be turned into positive ones.
 173
 174     3) Compute the even-mode impedance, assuming a vacuum dielectric, or
 175     if there is just one dielectric, that one.
 176
 177     4) Compute the even-mode impedance, taking into account the effect of
 178     multiple dielectrics, IF NECESSARY  */
 179
 180     /* Stage 1 - compute the odd mode impedance assuming single dielectric */
 181     data->display = Z_ODD_SINGLE_DIELECTRIC;
 182     dielectrics_to_consider_just_now=1;
 183     data->dielectrics_to_consider_just_now=1;
 184
 185     capacitance=VERY_LARGE; /* Can be anything large */
 186     if(data->verbose_level >= 2)
 187       printf("Solving assuming a vacuum dielectric to compute the odd-mode impedance\n");
 188
 189     do /* Start a finite difference calculation */
 190     {
 191       capacitance_old=capacitance;
 192 #ifdef ENABLE_POSIX_THREADS
 193       if (number_of_workers == 0)
 194          capacitance=finite_difference_single_threaded();
 195       else
 196         capacitance=finite_difference_multi_threaded();
 197 #else
 198       capacitance=finite_difference_single_threaded();
 199 #endif
 200       data->Codd_vacuum=capacitance;
 201       data->Codd=capacitance;
 202       data->Lodd_vacuum=MU_0*EPSILON_0/capacitance; /* Same as L in *ALL* cases */
 203
 204       data->Zodd_vacuum=sqrt(data->Lodd_vacuum/data->Codd_vacuum);  /* Standard formaul for Zodd */
 205
 206       if (data->dielectrics_in_bitmap == 1) /* Just get C by simple scaling of Er */
 207         data->Codd*=data->found_this_dielectric;  /* Scaled by the single dielectric constant */
 208       else
 209                 data->Er=1.0;
 210       data->Zodd=sqrt(data->Lodd_vacuum/data->Codd);  /* Standard formula for Zo */
 211       velocity_of_light_in_vacuum=1.0/(sqrt(MU_0 * EPSILON_0)); /* around 3x10^8 m/s */
 212       /* FPE trapdata->velocity_odd=1.0/pow(data->L_vacuum*data->Codd,0.5); */
 213       data->velocity_odd=1.0/pow(data->Lodd_vacuum*data->Codd,0.5);
 214       data->velocity_factor_odd=data->velocity_odd/velocity_of_light_in_vacuum;
 215       data->relative_permittivity_odd=sqrt(data->velocity_factor_odd); /* ??? XXXXXX */
 216       data->Er_odd=data->Codd/data->Codd_vacuum;
 217       data->Zdiff=2.0*data->Zodd;
 218       /* Print text if uses wants it */
 219       if(data->verbose_level>=1)
 220         print_data_for_directional_couplers(*data, where_to_print_fp, inputfile_filename);
 221     } while (fabs((capacitance_old-capacitance)/capacitance_old) > data->cutoff); /* end of FD loop */
 222
 223 #ifdef ENABLE_MPI
 224         mpi_receive_updated_vij_strips();
 225 #endif /* ENABLE_MPI */
 226
 227     /* display bitpamps/binary files if this is the last odd-mode computation */
 228     if((data->write_binary_field_imagesQ == TRUE || data->write_bitmap_field_imagesQ == TRUE) && data->dielectrics_in_bitmap==1 )
 229       write_fields_for_directional_couplers(inputfile_filename, *data, size, ODD);
 230
 231     /* Stage 2 - compute the odd-mode impedance taking into account other dielectrics IF NECESSARY */
 232
 233     if ( data->dielectrics_in_bitmap >1)
 234     {
 235       if(data->verbose_level >= 2)
 236         printf("Now taking into account the permittivities of the different dielectrics to compute Zodd.\n");
 237       data->display = Z_ODD_SINGLE_DIELECTRIC;
 238       capacitance=VERY_LARGE; /* Can be anything large */
 239
 240       dielectrics_to_consider_just_now=2;
 241       data->dielectrics_to_consider_just_now=2;
 242
 243 #ifdef ENABLE_MPI
 244           mpi_send_current_data();
 245 #endif /* ENABLE_MPI */
 246
 247       do /* Start a finite calculation */
 248       {
 249         capacitance_old=capacitance;
 250 #ifdef ENABLE_POSIX_THREADS
 251       if (number_of_workers == 0)
 252          capacitance=finite_difference_single_threaded();
 253       else
 254         capacitance=finite_difference_multi_threaded();
 255 #else
 256       capacitance=finite_difference_single_threaded();
 257 #endif
 258         data->Codd=capacitance;
 259         data->Zodd=sqrt(data->Lodd_vacuum/data->Codd);  /* Standard formula for Zo */
 260         velocity_of_light_in_vacuum=1.0/(sqrt(MU_0 * EPSILON_0)); /* around 3x10^8 m/s */
 261         data->velocity_odd=1.0/pow(data->L_vacuum*data->C,0.5);
 262         data->velocity_factor_odd=data->velocity/velocity_of_light_in_vacuum;
 263         data->relative_permittivity_odd=sqrt(data->velocity_factor); /* ??? XXXXXX */
 264                 data->Er_odd=data->Codd/data->Codd_vacuum;
 265                 data->Zdiff=2.0*data->Zodd;
 266                 if(data->verbose_level>=1)
 267           print_data_for_directional_couplers(*data, where_to_print_fp, inputfile_filename);
 268       } while (fabs((capacitance_old-capacitance)/capacitance_old) > data->cutoff); /* end of FD loop */
 269
 270 #ifdef ENABLE_MPI
 271           mpi_receive_updated_vij_strips();
 272 #endif /* ENABLE_MPI */
 273
 274       if((data->write_binary_field_imagesQ == TRUE || data->write_bitmap_field_imagesQ == TRUE) && data->dielectrics_in_bitmap!=1 )
 275         write_fields_for_directional_couplers(inputfile_filename, *data, size, ODD);
 276     } /* end of stage 2 for couplers */
 277
 278     /* Stage 3 - compute the even-mode impedance assuming single dielectric */
 279
 280     /* Since we want the even mode impedance now, we swap all the -1V
 281     metallic conductors for +1V */
 282
 283     swap_conductor_voltages();
 284
 285     data->display = Z_EVEN_SINGLE_DIELECTRIC;
 286     dielectrics_to_consider_just_now=1;
 287     data->dielectrics_to_consider_just_now=1;
 288     if(data->verbose_level >= 2)
 289         printf("Now assuming a vacuum dielectric to compute Zeven\n");
 290
 291     capacitance=VERY_LARGE; /* Can be anything large */
 292
 293 #ifdef ENABLE_MPI
 294         mpi_send_current_data();
 295 #endif /* ENABLE_MPI */
 296
 297     do /* Start a finite difference calculation */
 298     {
 299       capacitance_old=capacitance;
 300 #ifdef ENABLE_POSIX_THREADS
 301       if (number_of_workers == 0)
 302          capacitance=finite_difference_single_threaded();
 303       else
 304         capacitance=finite_difference_multi_threaded();
 305 #else
 306       capacitance=finite_difference_single_threaded();
 307 #endif
 308
 309       data->Ceven_vacuum=capacitance;
 310       data->Ceven=capacitance;
 311       data->Leven_vacuum=MU_0*EPSILON_0/capacitance; /* Same as L in *ALL* cases */
 312
 313       data->Zeven_vacuum=sqrt(data->Leven_vacuum/data->Ceven_vacuum);  /* Standard formaul for Zodd */
 314
 315       if (data->dielectrics_in_bitmap == 1) /* Just get C by simple scaling of Er */
 316         data->Ceven*=data->found_this_dielectric;  /* Scaled by the single dielectric constant */
 317       else
 318                 data->Er_even=1.0;
 319       data->Zeven=sqrt(data->Leven_vacuum/data->Ceven);  /* Standard formula for Zo */
 320       velocity_of_light_in_vacuum=1.0/(sqrt(MU_0 * EPSILON_0)); /* around 3x10^8 m/s */
 321       data->velocity_even=1.0/pow(data->Leven_vacuum*data->Ceven,0.5);
 322       data->velocity_factor_even=data->velocity_even/velocity_of_light_in_vacuum;
 323       data->relative_permittivity_even=sqrt(data->velocity_factor_even); /* ??? XXXXXX */
 324       data->Er_even=data->Ceven/data->Ceven_vacuum;
 325       data->Zcomm=data->Zeven/2.0;
 326       data->Zo=sqrt(data->Zodd * data->Zeven);
 327           if(data->verbose_level>=1)
 328                 print_data_for_directional_couplers(*data, where_to_print_fp, inputfile_filename);
 329       /* display bitpamps/binary files if this is the last even-mode computation */
 330     } while (fabs((capacitance_old-capacitance)/capacitance_old) > data->cutoff); /* end of FD loop */
 331
 332     if((data->write_binary_field_imagesQ == TRUE || data->write_bitmap_field_imagesQ == TRUE) && data->dielectrics_in_bitmap==1)
 333       write_fields_for_directional_couplers(inputfile_filename, *data, size, EVEN);
 334
 335     capacitance=VERY_LARGE; /* Can be anything large */
 336     /* Stage 4 - compute the even-mode impedance assuming multiple dielectrics IF NECESSARY */
 337     if ( data->dielectrics_in_bitmap >1)
 338     {
 339       dielectrics_to_consider_just_now=2;
 340       data->dielectrics_to_consider_just_now=2;
 341       if(data->verbose_level >= 2)
 342         printf("Now taking into account the permittivities of the different dielectrics to compute Zeven\n");
 343
 344 #ifdef ENABLE_MPI
 345           mpi_send_current_data();
 346 #endif /* ENABLE_MPI */
 347
 348       do /* Start a finite calculation */
 349       {
 350         capacitance_old=capacitance;
 351 #ifdef ENABLE_POSIX_THREADS
 352         if (number_of_workers == 0)
 353            capacitance=finite_difference_single_threaded();
 354         else
 355           capacitance=finite_difference_multi_threaded();
 356 #else
 357         capacitance=finite_difference_single_threaded();
 358 #endif
 359         data->Ceven=capacitance;
 360         data->Zeven=sqrt(data->Leven_vacuum/data->Ceven);  /* Standard formula for Zo */
 361         velocity_of_light_in_vacuum=1.0/(sqrt(MU_0 * EPSILON_0)); /* around 3x10^8 m/s */
 362         data->velocity_even=1.0/pow(data->L_vacuum*data->C,0.5);
 363         data->velocity_factor_even=data->velocity/velocity_of_light_in_vacuum;
 364         data->relative_permittivity_even=sqrt(data->velocity_factor); /* ??? XXXXXX */
 365                 data->Er_even=data->Ceven/data->Ceven_vacuum;
 366                 data->Zdiff=2.0*data->Zodd;
 367         data->Zcomm=data->Zeven/2.0;
 368                 data->Zo=sqrt(data->Zeven*data->Zodd);
 369                 if(data->verbose_level>=1)
 370           print_data_for_directional_couplers(*data, where_to_print_fp, inputfile_filename);
 371       } while (fabs((capacitance_old-capacitance)/capacitance_old) > data->cutoff); /* end of FD loop */
 372
 373       if(data->write_binary_field_imagesQ == TRUE || data->write_bitmap_field_imagesQ == TRUE)
 374         write_fields_for_directional_couplers(inputfile_filename, *data, size, EVEN);
 375     } /* end of stage 4 */
 376     /* Print the results if the verbose level was 0 (no -v flag(s) ). */
 377     if (data->verbose_level == 0)
 378     {
 379       /* We need to print the data. The next function will only print if
 380       the verbose_level is 1 or more, so I'll fix it at one. Then we print
 381       the final results and exit. */
 382       data->verbose_level=1;
 383       data->display = Z_EVEN_SINGLE_DIELECTRIC;
 384       print_data_for_directional_couplers(*data, where_to_print_fp, inputfile_filename);
 385     }
 386   } /* end of if couplers */
 387 }