net_solve.cpp

矩量法仿真电磁辐射和散射的源代码（c++）

						}
					}
			
					if ( ! jump2 )
					{
						irow1= ndimn- (ntsc+1);
						ntsca[ntsc]= nseg1;
						vsrc[ntsc]= source_voltage_array[isc1];
						ntsc++;
					}
			
				} /* if ( ! jump1 ) */
		
				jump1 = false;
				for( i = 0; i < voltage_source_count; i++ )
				{
					if ( nseg2 == source_segment_array[i])
					{
						isc2= i;
						jump1 = true;
						break;
					}
				}
		
				jump2 = false;
				if ( ! jump1 )
				{
					isc2=-1;
			
					for( i = 0; i < nteq; i++ )
					{
						if ( nseg2 == nteqa[i])
						{
							irow2= i;
							jump2 = true;
							break;
						}
					}
			
					if ( ! jump2 )
					{
						irow2= nteq;
						nteqa[nteq]= nseg2;
						nteq++;
					}
				}  /* if ( ! jump1 ) */
				else
				{
					for( i = 0; i < ntsc; i++ )
					{
						if ( nseg2 == ntsca[i])
						{
							irow2 = ndimn- (i+1);
							jump2 = true;
							break;
						}
					}
		
					if ( ! jump2 )
					{
						irow2= ndimn- (ntsc+1);
						ntsca[ntsc]= nseg2;
						vsrc[ntsc]= source_voltage_array[isc2];
						ntsc++;
					}
				} /* if ( ! jump1 ) */
		
				/* fill network equation matrix and right hand side vector with */
				/* network short-circuit admittance matrix coefficients. */
				if ( isc1 == -1)
				{
					cmn[irow1+irow1*ndimn] -= nec_complex( y11r, y11i)* geometry.segment_length[nseg1-1];
					cmn[irow1+irow2*ndimn] -= nec_complex( y12r, y12i)* geometry.segment_length[nseg1-1];
				}
				else
				{
					rhnx[irow1] += nec_complex( y11r, y11i)* source_voltage_array[isc1]/wavelength;
					rhnx[irow2] += nec_complex( y12r, y12i)* source_voltage_array[isc1]/wavelength;
				}
			
				if ( isc2 == -1)
				{
					cmn[irow2+irow2*ndimn] -= nec_complex( y22r, y22i)* geometry.segment_length[nseg2-1];
					cmn[irow2+irow1*ndimn] -= nec_complex( y12r, y12i)* geometry.segment_length[nseg2-1];
				}
				else
				{
					rhnx[irow1] += nec_complex( y12r, y12i)* source_voltage_array[isc2]/wavelength;
					rhnx[irow2] += nec_complex( y22r, y22i)* source_voltage_array[isc2]/wavelength;
				}
			} /* for( j = 0; j < network_count; j++ ) */
	
			/*	add interaction matrix admittance
				elements to network equation matrix */
			for( i = 0; i < nteq; i++ )
			{
				for (int j = 0; j < neqt; j++ )
					rhs[j] = cplx_00();
				
				irow1= nteqa[i]-1;
				rhs[irow1]=cplx_10();
				solves( cm, ip, rhs, neq, 1, geometry.np, geometry.n, geometry.mp, geometry.m, nop, symmetry_array);
				geometry.get_current_coefficients(wavelength, rhs, air, aii, bir, bii, cir, cii, vqds, nqds, iqds);
				
				for (int j = 0; j < nteq; j++ )
				{
					irow1= nteqa[j]-1;
					cmn[i+j*ndimn] += rhs[irow1];
				}
			} /* for( i = 0; i < nteq; i++ ) */
		
			/* factor network equation matrix */
			lu_decompose( nteq, cmn, ipnt, ndimn);
			
		} /* if ( network_count != 0) */
	} /* if ( ntsol != 0) */

	if (0 == network_count)
	{
		/* solve for currents when no networks are present */
		solves( cm, ip, einc, neq, 1, geometry.np, geometry.n, geometry.mp, geometry.m, nop, symmetry_array);
		geometry.get_current_coefficients(wavelength, einc, air, aii, bir, bii, cir, cii, vqds, nqds, iqds);
		ntsc=0;
	}
	else // if ( network_count != 0)
	{
		/* add to network equation right hand side */
		/* the terms due to element interactions */
		for( i = 0; i < neqt; i++ )
			rhs[i]= einc[i];
	
		solves( cm, ip, rhs, neq, 1, geometry.np, geometry.n, geometry.mp, geometry.m, nop, symmetry_array);
		geometry.get_current_coefficients(wavelength, rhs, air, aii, bir, bii, cir, cii, vqds, nqds, iqds);
	
		for( i = 0; i < nteq; i++ )
		{
			irow1= nteqa[i]-1;
			rhnt[i]= rhnx[i]+ rhs[irow1];
		}

		/* solve network equations */
		solve( nteq, cmn, ipnt, rhnt, ndimn);
	
		/* add fields due to network voltages to electric fields */
		/* applied to structure and solve for induced current */
		for( i = 0; i < nteq; i++ )
		{
			irow1= nteqa[i]-1;
			einc[irow1] -= rhnt[i];
		}
	
		solves( cm, ip, einc, neq, 1, geometry.np, geometry.n, geometry.mp, geometry.m, nop, symmetry_array);
		geometry.get_current_coefficients(wavelength, einc, air, aii, bir, bii, cir, cii, vqds, nqds, iqds);

		if ( nprint == 0)
		{
			fprintf( output_fp, "\n\n\n"
				"                          "
				"--------- STRUCTURE EXCITATION DATA AT NETWORK CONNECTION POINTS --------" );
		
			fprintf( output_fp, "\n"
				"  TAG   SEG       VOLTAGE (VOLTS)          CURRENT (AMPS)        "
				" IMPEDANCE (OHMS)       ADMITTANCE (MHOS)     POWER\n"
				"  No:   No:     REAL      IMAGINARY     REAL      IMAGINARY    "
				" REAL      IMAGINARY     REAL      IMAGINARY   (WATTS)" );
		}

		for( i = 0; i < nteq; i++ )
		{
			int segment_number = nteqa[i];
			int segment_index = segment_number-1;
			nec_complex voltage = rhnt[i]* geometry.segment_length[segment_index]* wavelength;
			nec_complex current = einc[segment_index]* wavelength;
			nec_complex admittance = current / voltage;
			nec_complex impedance = voltage / current;
			int segment_tag = geometry.segment_tags[irow1];
			nec_float power = em::power(voltage,current);
			network_power_loss= network_power_loss - power;
			
			if ( nprint == 0)
				fprintf( output_fp, "\n"
					" %4d %5d %11.4E %11.4E %11.4E %11.4E"
					" %11.4E %11.4E %11.4E %11.4E %11.4E",
					segment_tag, segment_number, real(voltage), imag(voltage), real(current), imag(current),
					real(impedance), imag(impedance), real(admittance), imag(admittance), power );
		}

		for( i = 0; i < ntsc; i++ )
		{
			irow1= ntsca[i]-1;
			vlt= vsrc[i];
			cux= einc[irow1]* wavelength;
			ymit= cux/ vlt;
			zped= vlt/ cux;
			irow2= geometry.segment_tags[irow1];
			
			nec_float pwr= em::power(vlt,cux);
			network_power_loss= network_power_loss- pwr;
		
			if ( nprint == 0)
				fprintf( output_fp, "\n"
					" %4d %5d %11.4E %11.4E %11.4E %11.4E"
					" %11.4E %11.4E %11.4E %11.4E %11.4E",
					irow2, irow1+1, real(vlt), imag(vlt), real(cux), imag(cux),
					real(zped), imag(zped), real(ymit), imag(ymit), pwr );
		} /* for( i = 0; i < ntsc; i++ ) */
	} /* if ( network_count != 0) */

	if ( (voltage_source_count+nvqd) == 0)
		return;
	
	nec_antenna_input* antenna_input = new nec_antenna_input();
	s_results.add(antenna_input);
	
	s_output.end_section();
	fprintf( output_fp, 
		"                        "
		"--------- ANTENNA INPUT PARAMETERS ---------" );
	
	fprintf( output_fp, "\n"
		"  TAG   SEG       VOLTAGE (VOLTS)         "
		"CURRENT (AMPS)         IMPEDANCE (OHMS)    "
		"    ADMITTANCE (MHOS)     POWER\n"
		"  NO.   NO.     REAL      IMAGINARY"
		"     REAL      IMAGINARY     REAL      "
		"IMAGINARY    REAL       IMAGINARY   (WATTS)" );
	
	for( i = 0; i < voltage_source_count; i++ )
	{
		int segment_index = source_segment_array[i]-1;
		nec_complex voltage = source_voltage_array[i];
		nec_complex current = einc[segment_index] * wavelength;
		
		bool add_as_network_loss = false;
		
		// the following loop is completely mysterious!
		for (int j = 0; j < ntsc; j++ )
		{
			// I am now almost sure that the following code is not correct.
			// This modifies the current, however if the inner loop is executed more
			// than once, then only the last current modification is kept!
			 
			if ( ntsca[j] == segment_index+1)
			{
				int row_index = ndimn - (j+1);
				int row_offset = row_index*ndimn;
				
				// I wish I knew what was going on here...
				nec_complex temp = rhnx[row_index]; // renamed current -> temp to avoid confusion
				for (int k = 0; k < nteq; k++ )
					temp -= cmn[k + row_offset]*rhnt[k];
					
				current = (temp + einc[segment_index])* wavelength;
				add_as_network_loss = true;
					
#warning "This loop is messed up. The j is inside another j loop"
				// I have removed the j from the "for (int k = 0; k < nteq; k++ )" loop 
				// and placed this"j=nteq" statement here.
				j = nteq;
			}
		}
			
		nec_complex admittance = current / voltage;
		nec_complex impedance = voltage / current;
		nec_float power = em::power(voltage,current);
		
		if ( add_as_network_loss )
			network_power_loss += power;
			
		input_power += power;
		
		int segment_tag = geometry.segment_tags[segment_index];
	
		antenna_input->set_input(
			segment_tag, segment_index+1,
			voltage, current, impedance, admittance, power);
		
		fprintf( output_fp,	"\n"
			" %4d %5d %11.4E %11.4E %11.4E %11.4E"
			" %11.4E %11.4E %11.4E %11.4E %11.4E",
			segment_tag, segment_index+1, real(voltage), imag(voltage), real(current), imag(current),
			real(impedance), imag(impedance), real(admittance), imag(admittance), power );
		
	} /* for( i = 0; i < voltage_source_count; i++ ) */

	
	for( i = 0; i < nvqd; i++ )
	{
		int segment_index = ivqd[i]-1;
		nec_complex voltage = vqd[i];
		
		nec_complex _ai( air[segment_index], aii[segment_index]);
		nec_complex _bi( bir[segment_index], bii[segment_index]);
		nec_complex _ci( cir[segment_index], cii[segment_index]);
		
		// segment length is measured in wavelengths. The pase is therefore the length in wavelengths
		// multiplied by pi().
		nec_float segment_length_phase = geometry.segment_length[segment_index] * pi(); // TCAM CHANGED TO pi() (from TP*.5)!!
		
		nec_complex current = ( _ai - _bi* sin(segment_length_phase)+ _ci * cos(segment_length_phase)) * wavelength;
		
		nec_complex admittance = current / voltage;
		nec_complex impedance = voltage / current;
		nec_float power = em::power(voltage,current);
		
		input_power += power;
		
		int segment_tag = geometry.segment_tags[segment_index];
	
		antenna_input->set_input(
			segment_tag, segment_index+1,
			voltage, current, impedance, admittance, power);
		
		fprintf( output_fp,	"\n"
			" %4d %5d %11.4E %11.4E %11.4E %11.4E"
			" %11.4E %11.4E %11.4E %11.4E %11.4E",
			segment_tag, segment_index+1, real(voltage), imag(voltage), real(current), imag(current),
			real(impedance), imag(impedance), real(admittance), imag(admittance), power );
		
	} /* for( i = 0; i < nvqd; i++ ) */
}