added 4nec2 LD 6 support (experimental!)

Author: maurymarkowitz

doc/TODO

@@ -1,7 +1,5 @@
 Items still to be implemented:
 
-* add support for the nec-2 style calling convention, using a tmp file with the input and output file names on separate lines
-
 * finish maa support, need to properly handle the loading based on wire and segment
 
 * get more info on the .ez format
@@ -22,7 +20,7 @@ Items still to be implemented:
 
 * nec2++ examples/medium_test.nec has an "MP" card type
 
-* at some point, handle 4nec2 current sources, EX 6 and 7 (UA3SFH.nec, W2EEY.NEC, calvin.nec, class.NEC)
+* at some point, handle 4nec2 current sources, EX 6 (UA3SFH.nec, W2EEY.NEC, calvin.nec, class.NEC)
 
 * add support for LD type 7 from 4nec2
 

src/calculations.c

@@ -351,9 +351,9 @@ int apply_impedance_loading(nec_context_t *ctx, int *ldtyp, int *ldtag, int *ldt
 
     } /* if( istep > ctx->zload.num_loads) */
 
-    /* ldtyp is validated (0–5) by control.c before storage; this should
+    /* ldtyp is validated (0–6) by control.c before storage; this should
      * never fire. If it does, it is a programming error, not a user error. */
-    assert(ldtyp[istepx] <= 5 && "INTERNAL: IMPROPER LOAD TYPE stored in zload.load_types");
+    assert(ldtyp[istepx] <= 6 && "INTERNAL: IMPROPER LOAD TYPE stored in zload.load_types");
 
     /* search segments for proper itags */
     ldtags= ldtag[istepx];
@@ -430,7 +430,30 @@ int apply_impedance_loading(nec_context_t *ctx, int *ldtyp, int *ldtag, int *ldt
 
         case 6:
           wire_surface_impedance( ctx, zlr[istepx]* ctx->geometry.wavelength, ctx->geometry.radius[i], &zt );
-          
+          break;
+
+        case 7: {
+          /* LD type 6: LC-trap.  Convert to parallel-RLC at the design frequency.
+           * F1 = unloaded-Q of the inductor (0 -> 100), F2 = L (H), F3 = C (F)
+           * F4 = optional design frequency in MHz (0 -> use first FR card). */
+          double f0_mhz = (ctx->zload.load_freq[istepx] != 0.0)
+                              ? ctx->zload.load_freq[istepx]
+                              : ctx->save.first_fr_mhz;
+          double Q     = (fabs(zlr[istepx]) < 1.0e-20) ? 100.0 : zlr[istepx];
+          double L     = zli[istepx];
+          double C     = zlc[istepx];
+          double omega0 = 2.0 * M_PI * f0_mhz * 1.0e6;
+          double R_par  = Q * omega0 * L;  /* parallel loss resistance at design freq */
+          /* Apply parallel-RLC formula (same as type 1 / case 2) */
+          zt = tpcj * ctx->geometry.half_len[i] * C / ctx->geometry.wavelength;
+          if (fabs(L) > 1.0e-20)
+            zt += ctx->geometry.half_len[i] * ctx->geometry.wavelength / (tpcj * L);
+          if (R_par > 1.0e-20)
+            zt += ctx->geometry.half_len[i] / R_par;
+          zt = 1.0 / zt;
+          break;
+        }
+
       } /* switch( jump ) */
 
       if(( fabs( creal( ctx->zload.seg_impedance[i]))+ fabs( cimag( ctx->zload.seg_impedance[i]))) > 1.0e-20) {
@@ -499,6 +522,10 @@ int apply_impedance_loading(nec_context_t *ctx, int *ldtyp, int *ldtag, int *ldt
 
       case 6:
            add_loading_output(ctx, ldtags, ldtagf[istepx], ldtagt[istepx], zlr[istepx], 0.0, 0.0, "WIRE");
+        break;
+
+      case 7:
+           add_loading_output(ctx, ldtags, ldtagf[istepx], ldtagt[istepx], 0.0, zli[istepx], zlc[istepx], "LC-TRAP");
 
     } /* switch( jump ) */
   } /* while( true ) */

src/control.c

@@ -434,6 +434,7 @@ static void reset_loading_buffers(nec_context_t *ctx)
         mem_free(ctx, (void **)&ctx->zload.load_r);
         mem_free(ctx, (void **)&ctx->zload.load_l);
         mem_free(ctx, (void **)&ctx->zload.load_c);
+        mem_free(ctx, (void **)&ctx->zload.load_freq);
         ctx->zload.num_loads = 0;
     }
     if (ctx->loading_outputs.entries != NULL) {
@@ -643,14 +644,16 @@ static int process_next_batch(nec_context_t *ctx, deck_t *deck, int *batch_start
             ctx->save.num_freq = (i2 == 0) ? 1 : i2;
             ctx->save.freq_mhz = f1;
             ctx->save.freq_step = f2;
+            if (ctx->save.first_fr_mhz == 0.0)
+                ctx->save.first_fr_mhz = f1;
         }
         else if (strcmp(code, "LD") == 0) {
             // LD card - Loading
             if (i1 == -1) {
                 continue;
             }
 
-            if (i1 > 5) {
+            if (i1 > 6) {
                 char msg[MAX_ERROR_LEN];
                 snprintf(msg, sizeof(msg), "LD on line %d: type %d is not supported.", card_idx + 1, i1);
                 add_error(ctx, &ctx->errors, msg, FATAL);
@@ -676,6 +679,7 @@ static int process_next_batch(nec_context_t *ctx, deck_t *deck, int *batch_start
             mem_realloc(ctx, (void **)&ctx->zload.load_r, mreq);
             mem_realloc(ctx, (void **)&ctx->zload.load_l, mreq);
             mem_realloc(ctx, (void **)&ctx->zload.load_c, mreq);
+            mem_realloc(ctx, (void **)&ctx->zload.load_freq, mreq);
 
             int idx = ctx->zload.num_loads - 1;
             ctx->zload.load_types[idx] = i1;
@@ -696,6 +700,9 @@ static int process_next_batch(nec_context_t *ctx, deck_t *deck, int *batch_start
             ctx->zload.load_r[idx] = f1;
             ctx->zload.load_l[idx] = f2;
             ctx->zload.load_c[idx] = f3;
+            /* For LD type 6: f4 optionally overrides the design frequency (MHz).  
+             * 0 means "use the first FR card's frequency" (resolved at compute time). */
+            ctx->zload.load_freq[idx] = (i1 == 6) ? f4 : 0.0;
         }
         else if (strcmp(code, "GN") == 0) {
             // GN card - Ground parameters  

src/deck_validations.c

@@ -561,7 +561,7 @@ void test_deck_structure(const nec_context_t *ctx, const deck_t *deck, errors_li
         int tag = deck->cards[i].i[2];
         int seg1 = deck->cards[i].i[3];
         int seg2 = deck->cards[i].i[4];
-        if (type < -1)
+        if (type < -1 || (type > 6 && type != -1))
         {
           snprintf(msg, sizeof(msg), "LD on line %d: unexpected type I1=%d.", i, type);
           add_error(ctx, errors, msg, 0);
@@ -588,6 +588,15 @@ void test_deck_structure(const nec_context_t *ctx, const deck_t *deck, errors_li
         snprintf(msg, sizeof(msg), "LD on line %d: load values F1..F3 are all zero.", i);
         add_error(ctx, errors, msg, 0);
       }
+      // LD type 6 (LC-trap): L and C are required
+      if (deck->cards[i].ints_used >= 1 && deck->cards[i].i[1] == 6)
+      {
+        if (deck->cards[i].flts_used < 3 || deck->cards[i].f[2] == 0.0 || deck->cards[i].f[3] == 0.0)
+        {
+          snprintf(msg, sizeof(msg), "LD on line %d (LC-trap): F2 (L) and F3 (C) must be non-zero.", i);
+          add_error(ctx, errors, msg, 0);
+        }
+      }
     }
     // record for open-end placement validation later (only for LD cards)
     if (strcmp(code, "LD") == 0 && deck->cards[i].ints_used >= 3)

src/internals.h

@@ -412,7 +412,8 @@ typedef struct
 		screen_wire_len,    /* scrwlt — Fortran SCRWLT: screen radial wire length (m) */
 		screen_wire_radius, /* scrwrt — Fortran SCRWRT: screen radial wire radius (m) */
 		freq_mhz,           /* fmhz — Fortran FMHZ: current frequency (MHz) */
-		freq_step;          /* delfrq: frequency step size */
+		freq_step,          /* delfrq: frequency step size */
+		first_fr_mhz;       /* frequency from the first FR card; used as LD6 design-freq default */
 
 } run_params_t;
 
@@ -588,9 +589,10 @@ typedef struct
 		*ldcard_num;	/* Deck line number of the originating LD card */
 
 	double
-		*load_r,  /* zlr — Fortran ZLR: R value (Ω, H, or F depending on type) */
-		*load_l,  /* zli — Fortran ZLI: L value */
-		*load_c;  /* zlc — Fortran ZLC: C value */
+		*load_r,    /* zlr — Fortran ZLR: R value (Ω, H, or F depending on type) */
+		*load_l,    /* zli — Fortran ZLI: L value */
+		*load_c,    /* zlc — Fortran ZLC: C value */
+		*load_freq; /* LD type-6 design frequency (MHz); 0 = use first FR card */
 
 	complex double *seg_impedance; /* zarray — Fortran ZARRAY: per-segment normalized impedance */
 } impedance_loading_t;

test/ld6_trap_test.deck

@@ -0,0 +1,12 @@
+CM LC-trap loading test (LD type 6)
+CM Based on example from 4nec2 manual: LD 6 2 1 16 16 5.5e-6 45e-12
+CM A simple dipole with an LC-trap on tag 2.
+CE
+GW 1 10 0 0 0 0.1 0 0 0.001
+GW 2 16 0.1 0 0 0.2 0 0 0.001
+GE 0
+FR 0 1 0 0 300
+LD 6 2 1 16 16 5.5e-6 45e-12
+EX 0 1 5 0 1 0
+RP 0 1 1 1000 90 0 0 0
+EN