/*  XNECVIEW - a program for visualizing NEC2 input and output data
 *
 *  Copyright (C) 1998-2006,2011,2020 Pieter-Tjerk de Boer -- pa3fwm@amsat.org
 *
 *  Distributed on the conditions of version 2 of the GPL: see the files
 *  README and COPYING, which accompany this source file.
 *
 *  This module parses the NEC input file (i.e., the file describing the
 *  antenna's structure).
 */

#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <ctype.h>

#include "xnecview.h"


/* some error codes for internal use: */
#define Err_misc          1
#define Err_scan          2
#define Err_nwires        3
#define Err_nsurfaces     4
#define Err_nexcis        5
#define Err_nloads        6
#define Err_nnetws        7


int nextabs;        /* next absolute segment number to be assigned */
char lastCard[2];   /* type of previous card */

double len3d(Point p1)    /* calculates length of 3D vector p1 */
{
   return sqrt(p1.x*p1.x+p1.y*p1.y+p1.z*p1.z);
}

Point int3d(Point p1, Point p2, double f)  /* calculate intermediate point (3D); f=0..1 -> p1..p2 */
{
   Point p;
   p.x=(1-f)*p1.x+f*p2.x;
   p.y=(1-f)*p1.y+f*p2.y;
   p.z=(1-f)*p1.z+f*p2.z;
   return p;
}

void out3d(Surface *su,double f)  /* calculates su->pa, with f determining the length of the arrow */
{
   Point d0,d1,no;
   d0.x=su->p0.x-su->pc.x;   d0.y=su->p0.y-su->pc.y;    d0.z=su->p0.z-su->pc.z;
   d1.x=su->p1.x-su->pc.x;   d1.y=su->p1.y-su->pc.y;    d1.z=su->p1.z-su->pc.z;
   no.x=d0.y*d1.z-d0.z*d1.y; no.y=-d0.x*d1.z+d0.z*d1.x; no.z=d0.x*d1.y-d0.y*d1.x;
   f *= (len3d(d0)+len3d(d1))/len3d(no);
   no.x*=f;                  no.y*=f;                   no.z*=f; 
   su->pa.x=su->pc.x+no.x;   su->pa.y=su->pc.y+no.y;    su->pa.z=su->pc.z+no.z;
}


void updateextremes(Point *p)
{
   if (p->x>extr_str) extr_str=p->x;
   if (p->y>extr_str) extr_str=p->y;
   if (p->z>extr_str) extr_str=p->z;
   if (-p->x>extr_str) extr_str=-p->x;
   if (-p->y>extr_str) extr_str=-p->y;
   if (-p->z>extr_str) extr_str=-p->z;
}


int reflect(int plane,int inc)           /* perform a reflection of all antenna elements */
{
   Wire *w,*w0;
   Surface *s,*s0;
   int i;

   EXPAND_IF_NECESSARY(numwires*2,maxwires,wires)
   for ( i=numwires, w0=wires, w=wires+numwires ; i>0 ; i--, w0++, w++, numwires++ ) {
      *w=*w0;
      switch (plane) {
         case 0: w->p0.x=-w->p0.x; w->p1.x=-w->p1.x; break;
         case 1: w->p0.y=-w->p0.y; w->p1.y=-w->p1.y; break;
         case 2: w->p0.z=-w->p0.z; w->p1.z=-w->p1.z; break;
      }
      if (w->itg>0) w->itg+=inc;
      w->abs=nextabs; nextabs+=w->ns;
   }
   EXPAND_IF_NECESSARY(numsurfaces*2,maxsurfaces,surfaces)
   for ( i=numsurfaces, s0=surfaces, s=surfaces+numsurfaces ; i>0 ; i--, s0++, s++, numsurfaces++ ) {
      *s=*s0;
      switch (plane) {
         case 0: s->p0.x=-s->p0.x; s->p1.x=-s->p1.x; s->p2.x=-s->p2.x; s->p3.x=-s->p3.x; s->pc.x=-s->pc.x; s->pa.x=-s->pa.x; break;
         case 1: s->p0.y=-s->p0.y; s->p1.y=-s->p1.y; s->p2.y=-s->p2.y; s->p3.y=-s->p3.y; s->pc.y=-s->pc.y; s->pa.y=-s->pa.y; break;
         case 2: s->p0.z=-s->p0.z; s->p1.z=-s->p1.z; s->p2.z=-s->p2.z; s->p3.z=-s->p3.z; s->pc.z=-s->pc.z; s->pa.z=-s->pa.z; break;
      }
   }
   return 0;
}


void removecommas(char *p)        /* change commas in string into spaces */
{
   while ((p=strchr(p,','))) *p++=' ';
}


int find_segment(int tag, int seg, char *card, Wire **wirep, int *segp)
{
   int i;
   if (tag>0) {                   /* location specified as tag,seg pair */
      for (i=0;i<numwires;i++) if (wires[i].itg==tag) break;
      if (i==numwires) { fprintf(stderr,"Unknown tag number on %s card: %i\n",card,tag); return Err_misc; }
      if (wires[i].ns>0) {
         if (seg<=0 || seg>wires[i].ns) { fprintf(stderr,"Incorrect segment number on %s card: %i\n",card,seg); return Err_misc; }
      } else {
         while (wires[i].ns!=-seg && wires[i].itg==tag) i++;
         if (wires[i].itg!=tag) { fprintf(stderr,"Incorrect segment number on %s card: %i\n",card,seg); return Err_misc; }
      }
      *segp = seg;
   } else {                       /* location specified as absolute segment number */
      for (i=0;i<numwires;i++) if (seg>=wires[i].abs && seg<wires[i].abs+wires[i].ns) break;
      if (i==numwires) { fprintf(stderr,"Incorrect absolute segment number on %s card: %i\n",card,seg); return Err_misc; }
      *segp = seg-wires[i].abs+1;
   }
   *wirep=wires+i;
   return 0;
}


int read_nec_CM(char *s)     /* CM -> comment card; check whether it contains xnecview options */
{
   s+=2;
   while (*s==' ') s++;
   if (memcmp(s,"xnecview:",9)==0) process_optionstring(s+9);
   return 0;
}


int read_nec_EX(char *s)     /* EX -> excitation */
{
   int type,tag,seg;
   int err;

   removecommas(s);
   if (sscanf(s,"EX%d%d%d",&type,&tag,&seg) != 3) return Err_scan;
   if (type!=0 && type!=5) { fprintf(stderr,"Only voltage sources can be displayed: %s\n",s); return Err_misc; }

   EXPAND_IF_NECESSARY(numexcis,maxexcis,excis)
   err=find_segment(tag, seg, "EX", &excis[numexcis].wire, &excis[numexcis].seg);
   if (err) return err;
   numexcis++;
   return 0;
}


int read_nec_TL(char *s)     /* TL -> transmission line */
{
   int tag1,seg1,tag2,seg2;
   int err;
   double imp;

   removecommas(s);
   if (sscanf(s,"TL%d%d%d%d%lg",&tag1,&seg1,&tag2,&seg2,&imp) != 5) return Err_scan;

   EXPAND_IF_NECESSARY(numnetws,maxnetws,netws)
   err=find_segment(tag1, seg1, "TL", &netws[numnetws].wire1, &netws[numnetws].seg1);
   if (err) return err;
   err=find_segment(tag2, seg2, "TL", &netws[numnetws].wire2, &netws[numnetws].seg2);
   if (err) return err;
   if (imp<0) netws[numnetws].type = -1; else netws[numnetws].type = 1;
   numnetws++;
   return 0;
}


int read_nec_NT(char *s)     /* NT -> network */
{
   int tag1,seg1,tag2,seg2;
   int err;

   removecommas(s);
   if (sscanf(s,"NT%d%d%d%d",&tag1,&seg1,&tag2,&seg2) != 4) return Err_scan;

   EXPAND_IF_NECESSARY(numnetws,maxnetws,netws)
   err=find_segment(tag1, seg1, "NT", &netws[numnetws].wire1, &netws[numnetws].seg1);
   if (err) return err;
   err=find_segment(tag2, seg2, "NT", &netws[numnetws].wire2, &netws[numnetws].seg2);
   if (err) return err;
   netws[numnetws].type = 0;
   numnetws++;
   return 0;
}


int addload(Wire *wi,int first,int last)   /* add one load */
{
   EXPAND_IF_NECESSARY(numloads,maxloads,loads)
   loads[numloads].wire=wi;
   loads[numloads].firstseg=first;
   loads[numloads].lastseg=last;
   numloads++;
   return 0;
}

int read_nec_LD(char *s)     /* LD -> loading */
{
   int type,tag,segf,segl;
   int i;
   Wire *wi;

   removecommas(s);
   tag=-2;
   segf=segl=0;
   if (sscanf(s,"LD%d%d%d%d",&type,&tag,&segf,&segl) < 1) return Err_scan;
   if (tag==-2) return Err_scan;

   if (type==-1) { numloads=0; return 0; }

   if (tag>0) {                   /* location specified as tag,seg pair */
      for (i=0;i<numwires;i++) if (wires[i].itg==tag) break;
      if (i==numwires) { fprintf(stderr,"Unknown tag number on LD card: %i\n",tag); return Err_misc; }
      if (wires[i].ns>0) {
         if (segf==0) segf=1;
         if (segl==0) segl=wires[i].ns;
         if (segf<=0 || segl>wires[i].ns || segf>segl) { fprintf(stderr,"Incorrect segment numbers on LD card: %i, %i\n",segf,segl); return Err_misc; }
         addload(wires+i,segf,segl);
         return 0;
      } else {
         if (segf==0) segf=-1;
         else segf=-segf;
         if (segl==0) segl=-INT_MAX;
         else segl=-segl;
         while (wires[i].ns!=segf && wires[i].itg==tag && i<numwires) i++;
         if (wires[i].itg!=tag || i>=numwires) { fprintf(stderr,"Incorrect first segment number on LD card: %i\n",segf); return Err_misc; }
         while (wires[i].ns>=segl && wires[i].itg==tag && i<numwires) {
            addload(wires+i,1,1);
            i++;
         }
         if ((wires[i].itg!=tag || i>=numwires) && segl!=-INT_MAX) { fprintf(stderr,"Incorrect last segment number on LD card: %i\n",-segl); return Err_misc; }
      }
   } else {                       /* location specified as absolute segment number */
      if (segf==0 && segl==0) {
         for (i=0;i<numwires;i++)  /* load all wires */
            if (wires[i].ns>0) addload(wires+i,1,wires[i].ns);
            else addload(wires+i,1,1);
         return 0;
      }
      if (segf==0 || segl==0 || segf>segl) { fprintf(stderr,"Incorrect absolute segment numbers on LD card: %i %i\n",segf,segl); return Err_misc; }
      for (wi=wires;wi<wires+numwires;wi++) {
         if (wi->ns>0) {
            int f,l;
            f=segf-wi->abs+1;
            if (f>wi->ns) continue;
            l=segl-wi->abs+1;
            if (l<1) continue;
            if (f<1) f=1;
            if (l>wi->ns) l=wi->ns;
            addload(wi,f,l);
         } else {
            if (wi->abs>=segf && wi->abs<=segl) addload(wi,1,1);
         }
      }
   }
   return 0;
}



int read_nec_GA(char *s)     /* GA -> wire arc */
{
   int ns,itg,segnr;
   double rada,ang1,ang2,rad;
   double phi,incphi;
   Point po;
   Wire *w;

   removecommas(s);
   if (sscanf(s,"GA%d%d%lg%lg%lg%lg",&itg,&ns,&rada,&ang1,&ang2,&rad) != 6) return Err_scan;
   EXPAND_IF_NECESSARY(numwires+ns,maxwires,wires)
   phi=ang1*M_PI/180;
   incphi=(ang2-ang1)*M_PI/180/ns;
   po.x=cos(phi)*rada; po.y=0; po.z=sin(phi)*rada; updateextremes(&po);
   w=wires+numwires; numwires+=ns;
   for (segnr=1; segnr<=ns; segnr++) {
      w->p0=po;
      phi+=incphi;
      po.x=cos(phi)*rada; po.y=0; po.z=sin(phi)*rada; updateextremes(&po);
      w->p1=po;
      w->rad=rad;
      w->itg=itg;
      w->ns=-segnr;
      w->abs=nextabs; nextabs++;
      w->dispnr=(segnr-1==ns/2);
      w++;
   }
   return 0;
}


int read_nec_GH(char *s)     /* GH -> helix */
{
   int ns,itg,segnr;
   double sp,hl,a1,b1,a2,b2,rad;
   double z,incz;
   double phi,incphi;
   double rx,incrx;
   double ry,incry;
   Point po;
   Wire *w;

   removecommas(s);
   if (sscanf(s,"GH%d%d%lg%lg%lg%lg%lg%lg%lg",&itg,&ns,&sp,&hl,&a1,&b1,&a2,&b2,&rad) != 9) return Err_scan;
   EXPAND_IF_NECESSARY(numwires+ns,maxwires,wires)
   if (hl==0) { fprintf(stderr,"Helix with length=0 not yet supported.\n"); return Err_misc; }
   z=0; incz=hl/ns;
   rx=a1; incrx=(a2-a1)/ns;
   ry=b1; incry=(b2-b1)/ns;
   phi=0;
   incphi=2*M_PI*incz/sp;
   if (hl<0) {
      incz=-incz;
      phi=M_PI/2;
   }
   po.x=cos(phi)*rx; po.y=sin(phi)*ry; po.z=z; updateextremes(&po);
   w=wires+numwires; numwires+=ns;
   for (segnr=1; segnr<=ns; segnr++) {
      w->p0=po;
      z+=incz; rx+=incrx; ry+=incry; phi+=incphi;
      po.x=cos(phi)*rx; po.y=sin(phi)*ry; po.z=z; updateextremes(&po);
      w->p1=po;
      w->rad=rad;
      w->itg=itg;
      w->ns=-segnr;
      w->abs=nextabs; nextabs++;
      w->dispnr=(segnr-1==ns/2);
      w++;
   }
   return 0;
}


void pointtransf(Point *p0,Point *p1,double *mat,double xs,double ys,double zs)
{
   Point po;
   po.x= xs + mat[0]*p0->x + mat[1]*p0->y + mat[2]*p0->z;
   po.y= ys + mat[3]*p0->x + mat[4]*p0->y + mat[5]*p0->z;
   po.z= zs + mat[6]*p0->x + mat[7]*p0->y + mat[8]*p0->z;
   *p1=po;
   updateextremes(&po);
}

int read_nec_GM(char *s)    /* GM -> coordinate transformation */
{
   int nrpt,its,i,itsi;
   double rox,roy,roz,xs,ys,zs;
   double mat[9],mmat[9];
   double co,si;
   Wire *w0,*w1,*wlast;
   Surface *s0,*s1;
   double its_double;

   removecommas(s);
   rox=roy=roz=xs=ys=zs=its=0;
   i=sscanf(s,"GM%d%d%lg%lg%lg%lg%lg%lg%lg",&itsi,&nrpt,&rox,&roy,&roz,&xs,&ys,&zs,&its_double);
   if (i<2) return Err_scan;

   its = (int)(its_double+0.5);   // round to double; the detour via double is needed because some programs write something like 1.1e1 instead of 11 in this field, and this is not forbidden according to the specification

   rox*=M_PI/180.0; roy*=M_PI/180.0; roz*=M_PI/180.0;

   mat[0]=mat[4]=mat[8]=1;
   mat[1]=mat[2]=mat[3]=0;
   mat[5]=mat[6]=mat[7]=0;

   co=cos(rox); si=sin(rox);
   memcpy(mmat,mat,sizeof(mat));
   mat[3] = co*mmat[3] - si*mmat[6];
   mat[4] = co*mmat[4] - si*mmat[7];
   mat[5] = co*mmat[5] - si*mmat[8];
   mat[6] = si*mmat[3] + co*mmat[6];
   mat[7] = si*mmat[4] + co*mmat[7];
   mat[8] = si*mmat[5] + co*mmat[8];

   co=cos(roy); si=-sin(roy);
   memcpy(mmat,mat,sizeof(mat));
   mat[0] = co*mmat[0] - si*mmat[6];
   mat[1] = co*mmat[1] - si*mmat[7];
   mat[2] = co*mmat[2] - si*mmat[8];
   mat[6] = si*mmat[0] + co*mmat[6];
   mat[7] = si*mmat[1] + co*mmat[7];
   mat[8] = si*mmat[2] + co*mmat[8];

   co=cos(roz); si=sin(roz);
   memcpy(mmat,mat,sizeof(mat));
   mat[0] = co*mmat[0] - si*mmat[3];
   mat[1] = co*mmat[1] - si*mmat[4];
   mat[2] = co*mmat[2] - si*mmat[5];
   mat[3] = si*mmat[0] + co*mmat[3];
   mat[4] = si*mmat[1] + co*mmat[4];
   mat[5] = si*mmat[2] + co*mmat[5];

   EXPAND_IF_NECESSARY(numwires*(nrpt+1),maxwires,wires)
   if (nrpt==0) {
      w1=wires;
   } else {
      w1=wires+numwires;
   }

   w0=wires;
   i=0;
   do {
      wlast=wires+numwires;

      while (w0<wlast) {
         if (w0->itg>=its) {
            pointtransf(&w0->p0,&w1->p0,mat,xs,ys,zs);
            pointtransf(&w0->p1,&w1->p1,mat,xs,ys,zs);
            w1->rad=w0->rad;
            w1->ns=w0->ns;
            if (w0->itg>0) w1->itg=w0->itg+itsi;
            else w1->itg=w0->itg;
            w1->abs=nextabs; nextabs+=w1->ns;
            w1->dispnr=w0->dispnr;
            w1++;
            numwires++;
         }
         w0++;
      }
   } while (++i < nrpt);

   if (nrpt==0) {
      s0=s1=surfaces;
      i=numsurfaces;
   } else {
      int oldnumsurfaces=numsurfaces;
      i=numsurfaces*nrpt;
      numsurfaces*=(1+nrpt);
      EXPAND_IF_NECESSARY(numsurfaces,maxsurfaces,surfaces)
      s0=surfaces;
      s1=surfaces+oldnumsurfaces;
   }
   while (i--) {
      s1->type=s0->type;
      pointtransf(&s0->p0,&s1->p0,mat,xs,ys,zs);
      pointtransf(&s0->p1,&s1->p1,mat,xs,ys,zs);
      pointtransf(&s0->p2,&s1->p2,mat,xs,ys,zs);
      pointtransf(&s0->p3,&s1->p3,mat,xs,ys,zs);
      pointtransf(&s0->pc,&s1->pc,mat,xs,ys,zs);
      pointtransf(&s0->pa,&s1->pa,mat,xs,ys,zs);
      s0++; s1++;
   }

   return 0;
}

int read_nec_GR(char *s)    /* GR -> generate cylindrical structure */
{
   char gm[60];
   int n,inc;
   removecommas(s);
   if (sscanf(s,"GR%d%d",&inc,&n) != 2) return Err_scan;
   sprintf(gm,"GM %i %i 0 0 %g 0 0 0",inc,n-1,360./n);
   read_nec_GM(gm);           /* GR card is translated into equivalent GM card */
   return 0;
}

int read_nec_GS(char *s)      /* GS -> geometry scale */
{
   double a;
   Wire *w;
   Surface *su;
   int i;
   removecommas(s);
   if (sscanf(s,"GS%*d%*d%lg",&a) != 1) return Err_scan;
   for (i=0,w=wires;i<numwires;i++,w++) {
      w->p0.x*=a;
      w->p0.y*=a;
      w->p0.z*=a;
      w->p1.x*=a;
      w->p1.y*=a;
      w->p1.z*=a;
      w->rad*=a;
   }
   for (i=0,su=surfaces;i<numsurfaces;i++,su++) {
      su->p0.x*=a; su->p0.y*=a; su->p0.z*=a;
      su->p1.x*=a; su->p1.y*=a; su->p1.z*=a;
      su->p2.x*=a; su->p2.y*=a; su->p2.z*=a;
      su->p3.x*=a; su->p3.y*=a; su->p3.z*=a;
      su->pc.x*=a; su->pc.y*=a; su->pc.z*=a;
      su->pa.x*=a; su->pa.y*=a; su->pa.z*=a;
   }
   extr_str*=a;
   return 0;                               
}

int read_nec_GW(char *s)    /* GW -> straight wire */
{
   Wire *w;
   int i;
   EXPAND_IF_NECESSARY(numwires,maxwires,wires)
   w=wires+numwires++;
   removecommas(s);
   w->rad=0;
   i=sscanf(s,"GW%d%d%g%g%g%g%g%g%lg",&w->itg,&w->ns,&w->p0.x,&w->p0.y,&w->p0.z,&w->p1.x,&w->p1.y,&w->p1.z,&w->rad);
   if (i!=8 && i!=9) return Err_scan;  /* note: omitting the radius seems to be allowed, for tapered wires */
   w->abs=nextabs; nextabs+=w->ns;
   w->dispnr=1;
   updateextremes(&w->p0);
   updateextremes(&w->p1);
   return 0;
}

int read_nec_GX(char *s)     /* GX -> reflection in coordinate planes */
{
   int i,inc, err;
   removecommas(s);
   if (sscanf(s,"GX%d%d",&inc,&i) != 2) return Err_scan;
   if (i%10==1)      { err=reflect(2,inc); if (err) return err; inc*=2; }
   if ((i/10)%10==1) { err=reflect(1,inc); if (err) return err; inc*=2; }
   if ((i/100)==1) return reflect(0,inc);
   return 0;
}

int read_nec_SC(char *s)     /* SC -> continuation of SM, SP, or SC */
{
   Surface *su;
   int ns;
   Point p3;
   int n;

   if(lastCard[0]!='S'){ fprintf(stderr,"SC card not preceded by SM, SP, or SC.\n"); return Err_misc; }

   EXPAND_IF_NECESSARY(numsurfaces,maxsurfaces,surfaces)
   su=surfaces+numsurfaces++;
   removecommas(s);

   switch (lastCard[1]) {
      case 'M':
         if (sscanf(s,"SC%*i%*i%g%g%g",&su->p2.x,&su->p2.y,&su->p2.z) != 3) return Err_scan;
         updateextremes(&su->p2);
         su->pc=int3d(su->p0,su->p2,0.5);
         p3=int3d(su->p1,su->pc,2);
         updateextremes(&p3);
         out3d(su,0.05);
         su->type=SU_rect;
         break;
      case 'P':
         n=sscanf(s,"SC%*i%*i%g%g%g%g%g%g",&su->p2.x,&su->p2.y,&su->p2.z,&su->p3.x,&su->p3.y,&su->p3.z);
         if (n!=6 && n!=3) return Err_scan;
         updateextremes(&su->p2);
         switch (su->type) {
            case SU_rect:
               su->pc=int3d(su->p0,su->p2,0.5);
               su->p3=int3d(su->p1,su->pc,2);    // we don't really need su->p3 now, but do need it in case of continuation using more SC cards
               updateextremes(&su->p3);
               out3d(su,0.05);
               break;
            case SU_tri:
               su->pc=int3d(int3d(su->p0,su->p1,0.5),su->p2,1./3);
               out3d(su,0.05);
               break;
            case SU_quad:
               if (n!=6) return Err_scan;
               updateextremes(&su->p3);
               su->pc=int3d(int3d(su->p0,su->p1,0.5),int3d(su->p2,su->p3,0.5),0.5);
               out3d(su,0.05);
               break;
            default:
               fprintf(stderr,"SC card for surface type which is not rectangle/triangle/quadrilateral.\n"); return Err_misc; 
         }
         break;
      case 'C':
         n=sscanf(s,"SC%*i%d%g%g%g%g%g%g",&ns,&su->p2.x,&su->p2.y,&su->p2.z,&su->p3.x,&su->p3.y,&su->p3.z);
         if (n!=7 && n!=4) return Err_scan;
         Surface *suLast=surfaces+(numsurfaces-2);
         su->p0=suLast->p3;
         su->p1=suLast->p2;
         updateextremes(&su->p2);
         switch (ns) {
            case 1:
               su->type=SU_rect;
               su->pc=int3d(su->p0,su->p2,0.5);
               su->p3=int3d(su->p1,su->pc,2);
               updateextremes(&su->p3);
               out3d(su,0.05);
               break;
            case 3:
               if (n!=7) return Err_scan;
               su->type=SU_quad;
               updateextremes(&su->p3);
               su->pc=int3d(int3d(su->p0,su->p1,0.5),int3d(su->p2,su->p3,0.5),0.5);
               out3d(su,0.05);
               break;
            default:
               fprintf(stderr,"SC card specifying unsupported shape NS=%d.\n",ns); return Err_misc; 
         }
         break;
      default:
         fprintf(stderr,"SC card following card %c%c.\n",lastCard[0],lastCard[1]); return Err_misc; 
   }

   return 0;
}

int read_nec_SM(char *s)     /* SM -> multiple surface patch (shown as one large rectangular piece of surface) */
{
   Surface *su;
   EXPAND_IF_NECESSARY(numsurfaces,maxsurfaces,surfaces)
   su=surfaces+numsurfaces;
   removecommas(s);
   if (sscanf(s,"SM%*i%*i%g%g%g%g%g%g",&su->p0.x,&su->p0.y,&su->p0.z,&su->p1.x,&su->p1.y,&su->p1.z) !=6) return Err_scan;
   updateextremes(&su->p0);
   updateextremes(&su->p1);
   return 0;
}

int read_nec_SP(char *s)     /* SP -> surface patch */
{
   Surface *su;
   int ns;
   double x1,y1,z1,x2,y2,z2;
   double r;

   EXPAND_IF_NECESSARY(numsurfaces,maxsurfaces,surfaces)
   su=surfaces+numsurfaces;  
   su->type=SU_arb;
   removecommas(s);
   if (sscanf(s,"SP%*d%d%lg%lg%lg%lg%lg%lg",&ns,&x1,&y1,&z1,&x2,&y2,&z2) !=7) return Err_scan;

   if (ns==0) {
      su->type=SU_arb;
      su->pc.x=x1; su->pc.y=y1; su->pc.z=z1;
      updateextremes(&su->pc);
      r=sqrt(z2/M_PI);
      x2*=M_PI/180; y2*=M_PI/180;
      su->p0.x = x1 - r*sin(y2);
      su->p0.y = y1 + r*cos(y2);
      su->p0.z = z1;
      su->p1.x = x1 - r*sin(x2)*cos(y2);
      su->p1.y = y1 - r*sin(x2)*sin(y2);
      su->p1.z = z1 + r*cos(x2);
      su->p2=int3d(su->p0,su->pc,2);
      out3d(su,0.1);
      ++numsurfaces;
      return 0;      
   }

   switch (ns) {
      case 1:
         su->type=SU_rect;
         break;
      case 2:
         su->type=SU_tri;
         break;
      case 3:
         su->type=SU_quad;
         break;
      default:
         fprintf(stderr,"Unknown patch shape: NS=%i.\n",ns); return Err_misc;
   }
   su->p0.x=x1; su->p0.y=y1; su->p0.z=z1;
   su->p1.x=x2; su->p1.y=y2; su->p1.z=z2;
   updateextremes(&su->p0);
   updateextremes(&su->p1);
   return 0;
}

void read_nec_input(FILE *f)               /* read the input file, and update the arrays of wires and surfaces accordingly */
{
   char s[100];
   int err;
   int last_was_NT_or_TL=0;

   nextabs=1;
   fgets(s,100,f);
   while (!feof(f) && (s[0]!='E' || s[1]!='N')) {
      s[0]=toupper(s[0]);
      s[1]=toupper(s[1]);
      err=0;
      switch (s[0]) {
         case 'C':
            if (s[1]=='M') err=read_nec_CM(s);
            break;
         case 'E': 
            if (s[1]=='X') err=read_nec_EX(s);
            break;
         case 'G': 
            switch (s[1]) {
               case 'A': err=read_nec_GA(s); break;
               case 'C': break;  /* just ignore GC cards: they only modify the radius of GW cards, which xnecview doesn't use anyway */
               case 'H': err=read_nec_GH(s); break;
               case 'M': err=read_nec_GM(s); break;
               case 'R': err=read_nec_GR(s); break;
               case 'W': err=read_nec_GW(s); break;
               case 'X': err=read_nec_GX(s); break;
               case 'S': err=read_nec_GS(s); break;
               case 'F':
                  fprintf(stderr,"Warning: unhandled cardtype in %s",s);
                  break;
            }
            break;
         case 'L': 
            if (s[1]=='D') err=read_nec_LD(s);
            break;
         case 'M': 
            switch (s[1]) {
               case 'M':
               case 'E':
                  printf("%s",s);
                  break;
            }
            break;
         case 'N': 
            if (s[1]=='T') {
               if (!last_was_NT_or_TL) numnetws=0;
               err=read_nec_NT(s);
               last_was_NT_or_TL=1;
            }
            break;
         case 'S':
            switch (s[1]) {
               case 'C': err=read_nec_SC(s); break;
               case 'M': err=read_nec_SM(s); break;
               case 'P': err=read_nec_SP(s); break;
            }
            break;
         case 'T': 
            if (s[1]=='L') {
               if (!last_was_NT_or_TL) numnetws=0;
               err=read_nec_TL(s);
               last_was_NT_or_TL=1;
            }
            break;
      }
      switch (err) {
         case Err_scan: fprintf(stderr,"Error in parameters: %s\n",s); break;
         case Err_nwires: fprintf(stderr,"Too many wires.\n"); break;
         case Err_nsurfaces: fprintf(stderr,"Too many surfaces.\n"); break;
         case Err_nexcis: fprintf(stderr,"Too many excitations.\n"); break;
         case Err_nloads: fprintf(stderr,"Too many loads.\n"); break;
         case Err_nnetws: fprintf(stderr,"Too many networks and/or transmission lines.\n"); break;
      }
      lastCard[0]=s[0];
      lastCard[1]=s[1];
      fgets(s,100,f);
   }
   extr = extr_str;
}